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The ground states and excitations of two-dimensional insulating and doped Mott insulators are 
described by a bond operator formalism. While the method represents the degrees of freedom of an 
arbitrary antiferromagnet exactly, it is especially suited to systems in which there is a natural pairing 
of sites into bonds, as in states with spontaneous or explicit spin-Peierls order (or bond-centered 
charge order). In the undoped insulator, as discussed previously, we obtain both paramagnetic and 
magnetically-ordered states. We describe the evolution of superconducting order in the ground state 
with increasing doping — at low doping, the superconductivity is weak, can co-exist with magnetic 
order, and there are no gapless spin 1/2 fermionic excitations; at high doping, the magnetic order 
is absent and we obtain a BCS d-wave superconductor with gapless spin 1/2, nodal fermions. We 
present the critical theory describing the onset of these nodal fermionic excitations. We discuss 
the evolution of the spin spectrum, and obtain regimes where a spin 1 exciton contributes a sharp 
resonance in the dynamic spin susceptiblity. We also discuss the experimental consequences of 
low-energy, dynamically fluctuating, spin-Peierls order in an isotropic Cu02 plane — we compute 
consequences for the damping and dispersion of an optical phonon involving primarily the O ions, 
and compare the results with recent neutron scattering measurements of phonon spectra. 



I. INTRODUCTION 

By now, it is reasonably well established that the 
doped antiferromagnets found in the cuprate compounds 
have a superconducting ground state with a d-wave pair- 
ing symmetry. Moreover, many low temperature (T) 
properties appear to be well described in the framework 
of the conventional BCS theory of c?-wave superconduc- 
tors. There are a number of fascinating properties at 
temperatures above Tc (the critical temperature for the 
onset of superconductivity) which are not well under- 
stood, but there are numerous plausible candidate theo- 
ries for these, involving crossovers between different com- 
peting orders in doped antiferromagnetad. 

However, there are some low T properties of the su- 
perconducting state that do not appear naturally in the 
traditional BCS framework. Among these are (a) the ap- 
pearance of a 5 = 1/2 moment near Mji-magnetic Zn or 
Li impurties in the underdoped regiona'Q, (6) the presence 
of low energy collective spin excitations (a S" = 1 spin ex- 
citon) at (tt, tt) and nearby incommensurate wavevectors, 
and (c) instabilities to various co-existing spin and charge 
density wave states. While it is possible to "cook-up" a 
microscopic Hamiltonian, and a corresponding Hartree- 
Fock treatment, to generate any of these physical proper- 
ties in the superconducting state, a proper understanding 
of the physics should require that they emerge naturally 
from some deeper principle. 

In a series of papersat2l, it has been proposed that all 
of these unusual properties of the superconductor emerge 
very naturally in a theory of the doped Mott insulator. 
The theory is best understood in two distinct steps. In 



the first step, we disrupt the magnetic Neel long-range or- 
der in the insulator by adding some frustrating exchange 
interactions. Then, in step two, we dope the paramag- 
netic Mott insulator obtained in step one with mobile 
charge carriers. 

Considerable theoretical progress has been made in un- 
derstanding the first step: the quantum transition in- 
volving the destruction of Neel long-range order in the 
half-filled insulator. It has been argued that the first 
paramagnetic state on the other side of such a quantum 
critical point generically has the following properties: 

(i) a stable S = \ particle excitation (an exciton) and 
permanent confinement of S* = 1/2 spinous, 

(ii) broken translational symmetry due to the appearance 
of a bond-centered charge {e.g. spin-Peierls) order, and 

(iii) confinement of a S' = 1/2 moment in the vicinity of 
non-magnetic Zn/Li impurities. 

Note the similarity of the features (i)-(iii) of the in- 
sulator to the properties (o)-(c) of the superconductor. 
The essence of step two is then very simple: properties 
(i)-(iii) survive for a finite range of doping the paramag- 
netic Mott insulator, and this provides a "natural" un- 
derstanding of properties (a)-(c) of the superconductor. 

We now present a somewhat more detailed discussion 
of step two: doping the confined Mott insulator with 
mobile charge carriers. Apart from a possible insulating 
Wigner crystal ground states at very low doping (such 
a Wigner crystal state must appear in the presence of 
long-range Coulomb interactional—the ground state is ex- 
pected to be superconductingJ'EjQ, and this will also be 
the case in the calculations in the present naper. The 
pairing amplitude is found to be d-wave likeErLl i.e. the 
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pairing amplitude has opposite signs in the horizontal 
and vertical directions. By adjusting the strength of the 
Coulomb repulsion between the holes, we can modify the 
characteristic size of the hole pairs. Only for the larger 
hole pairs does the superconducting ground state posses 
gapless nodal fermionic excitations. We discuss the fate 
of (i)-(iii) in the doped antiferromagnct, along with con- 
nections to the experiments in turn: 

(i) S = 1 collective spin exciton: In the confined insula- 
tor, the lowest energy excitation above the spin gap is 
a stable, = 1, bosonic exciton with a minimum in its 
dispersion at (tTjTt). Upon doping to the superconduc- 
tor, additional gapped-jS* = 1/2 fermionic excitations are 
expected to appeait3ll3 (these are the Bogoliubov quasi- 
particles). For sufficiently large doping, the fermionic 
spectral becomes gapless at four nodal points in the Bril- 
louin zone, and so the global spin gap vanishes. In gen- 
eral, the S — 1 boson will be unstable to decay in into 
two of these fermionic excitations. However, constraints 
from momentum conservation can (and do) protect the 
integrity of the 5=1 exciton in certain regions of the 
Brillouin zone. In particular, the S = 1 exciton is likely 
to be stable near momenta (tt, tt); this can happen even if 
the fermionic excitations are gapless, provided the spac- 
ing between the nodal points does not equal (tt, tt). Ex- 
perimentally a S' = 1 neutron scattering resonajjce-is in- 
deed observed in the d-wave superconductorstirlla, and 
it is our contention that this excitation is continuously 
connected to that in the spin-Peierls insulator. We will 
obtain explicit results for the evolution of the bosonic 
and fermionic spin excitations in this paper, from the 
spin-Peierls insulator, to the superconductor with gap- 
less nodal excitations. 

(ii) Broken translational symmetry: In the simplest 



scenaricOlJ, which can be realized for a range of couplings 
in models without long-range Coulomb interactions, the 
bond-centered charge (or spin Peierls) order of the para- 
magnetic insulator survives in the doped superconductor 
all the way upto a critical doping at which full square lat- 
tice symmetry is restored in a transition to the c?-wave su- 
perconductor; with Coulomb interactions, Wigner crystal 
states also appear for a range of very small hole doping, 
but otherwise the situation is similar. The critical theory 
of the vanishing of the bond-centered (or site-centered) 
CDW order in tlae d-wave superconductor has been dis- 
cussed elsewhereQ. For other parameters, more complex 
striped states are also possible, with a period larger than 
2 sites and with modulation of the hole density on the 
sites. Experimentally, chapgp-stripe states of period 4 
have been clearly observecill3c3, but it is not yet estab- 
lished whether the modulation of the spin density is site 
or bond centepd. More recently, the observation of Mc- 
Queeney et alB^ in optimally doped and superconducting 
LSCO have been interpreted using a picture of bond- 
centered charge stripes of period 2, like those found in 
the spin-Peierls state . We will discuss these observa- 
tions further in Section |^. 



(iii) 3=1/2 moment near Zn/Li impurities: In princi- 
ple, it is possible that the S ~ 1/2 moment confined 
near a Zn/Li impurity in the paramagnetic Mott in- 
sulator disappears immediately at an infinitesimal hole 
doping concentration: one hole can be trapped near 
the Zn/Li impurity, and this configuration is compati- 
ble with the global preservation of the spin-Pcicrls (or 
other bond-centered charge-) order. However the kinetic 
energy cost, makes this unlikely. Barring this uninter- 
esting possibility, the S = 1/2 moment will survive in 
the superconductor, and this offers a natural explanation 
for the NMR experimentsB. Eventually, the fermionic 
S" = 1/2 excitations. p£ the superconductor will Kondo 
screen the momcntES^O, but because of the linearly van- 
ishing fermionic density of states at the Fermi level, this 
happens only above finite values of the impurity exchange 
coupling and particle-hole asymmetry (i.e. above a crit- 
ical doping). Moreover, there is no fundamental reason 
for this Kondo screening transition to co-incide with the 
point at which translational symmetry is restored (this 
transition was discussed above in (i)); the two transitions 
could occur in either order as a function of increasing 
doping 

The purpose of this paper is to present a theory 
of doped antiferromagents which displays the crossover 
from an insulator at zero doping (with or without long- 
range magnetic order) to a superconductor with gapless, 
nodal, fermionic 5 = 1/2 excitations at some moderate 
doping. Further, we require the theory to obey prop- 
erties (i)-(iii) directly at the mean-field level. While a 
large number of previous theories of doped antiferromag- 
nets have been presented previously, none of them satisfy 
all of these requirements. The studies of Ref |^ were able 
to examine the intricate competition between different 
charge ordered states and superconductivity — however, 
the ground states were well away from a region of mag- 
netic order, and there was no sharp, collective, 5=1 
excitation in the Gaussian fluctuations about the mean- 
field theory. Conversely, approaches which do yield con- 
finement of spinons. and collective 5=1 excitations at 
zero or low dopingQ'EStJ are not easily extended to reach 
a superconducting state with 5=1/2 gapless nodal ex- 
citations at moderate doping. 

We now outline the remainder of the paper. In Sec- 
tion H we introduce the central formalism of bond opera- 
tors, and its application to doped antiferromagnets. The 
main results of bond operator theory of the square lattice 
antiferromagnct are presented in Section HI, with most 
details of the mean-field calculations being relegated to 
Appendices refsquare, and ^: the important phase 
diagram is in Fig ^. The critical theory for the onset 
of nodal fermion excitations is in Section |^ Section 
differs from the remainder of the paper in that it consid- 
ers systems in which the bond-centered charge order is 
not present. Instead it considers the case when the spin- 
Peierls order is dynamically fluctuating and describes its 
influence on the optical phonon spectrum of the Cu02 
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plane — we find that the calculated phonon spectrm is in- 
deed strikingly similar to recent neutron scattering obser- 
vations. The main conclusions are stated in Section [v| . 
Finally, in Appendix ^ we describe the interplay between 
the spin-Peierls states considered in this paper, and states 
with electronic "nematic" order. 



II. BOND OPERATORS 

Our approach is a generalization of the bond opera- 
tor theory of Ref ^ (and the related work a£ Ref |2|) 
to doped antiferromagnets. The earlier workEj was de- 
signed for insulating systems, and has since been ap- 
plied succesfuUy in a ruitr^r of studies of spin ladder 
and related compjCumdsEa E3, and also in bilayer quan- 
tum Hall system£3H . Here we shall extend the for- 
malism to doped antiferromagneta;-a closely related ex- 
tension was discussed by Lee eL-aZ.LJ but only applied to 
one-dimensional systems. Ederuj, SushkovliJ, and Vojta 
and Beckero have also considered doped systemsJpy as- 
sociated methods. Recently, Jurecka and Breniglj have 
used a bond operator method to analyze a Kondo lattice 
model using a formulation that bears some similarity to 
ours. 

While our bond operator formalism can, in principle, 
be applied to an arbitrary doped antiferromagnet, it in- 
duces a bias by using a basis of states which explicitly 
refer to a preferred, disjoint pairing of all the sites. A 
full and exact solution of the bond operator Hamiltonian 
should restore the full symmetry of the underlying Hamil- 
tonian in which this preferred pairing may be absent. 
However, in practice, this restoration of symmetry is dif- 
ficult to achieve, and this is the principal shortcoming of 
the bond operator method. So the main utility of the ap- 
proach lies in treating systems in which there is a natural 
pairing of sites in the ground state, either imposed by the 
Hamiltonian, or by a spontaneous symmetry breaking. 
In this paper, we will restrict consideration to the case 
where the ground state possesses bond-centered charge 
density wave (spin-Peierls) order: there is evidently a 
natural pairing of sites in such a structure. The spin- 
Peierls order can either be spontaneous or explicit; the 
latter is the case in the doped two-leg ladder compounds 
like SrCu203 and Sri4_2;Caa;Cu2404i. In the presence 
of this background spin Peierls order, the method can 
then address the competition between magnetic Neel or- 
der and superconductivity, and follow the evolution of the 
fermionic 5=1/2 and the bosonic S* = 1 excitations. 

We will now introduce the formalism by showing the 
exact mappings between bond and site operators of a pair 
of sites. In hole-doped antiferromagnets, we can project 
out all states with 2 electrons on one site, and this leaves 
a total of 9 states in a pair of sites. Let and 

=1, 1) be the electron creatipn operators on the two 
sites. Then, as in the insulatorEj, we introduce four bond 



boson creation operators, and ij, (a = x,y,z) which 



are defined by (cr^f, are the Pauli matrices, and Eab is the 
second-rank antisymmetric tensor with £|| — 1) 



4l«> = -^<^bc£cac\aclt\0), 



(1) 



where |0) is the electron vacuum i.e. the state with no 
electrons on the two sites, while \v) is an unphysical state 
in which none of the bond bosons or fermions are present. 
To describe the remaining 5 states of the doped antiferro- 
magnet, we introduce the bond fermionic operators 
and /ija, and the additional bond bosonic operator 
which are defined by 



(2) 



The operators s, d, all obey the canonical boson com- 
mutation relations, while the hia, ft.2a obey canonical 
fermion relations. Of course, the total space of states 
in Fock space of these 5 bosons and 4 fermions is much 
larger than the 9 states allowed in the doped antiferro- 
magnet. To restrict to the physical subspace we must 
impose the single constraint 
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= cL|0) 






d^\v) 


= |0). 



S^S + tlta + h\jlia + hl^h2a + dU = 1. 



(3) 



In the subspace constrained by (|^), we can now write 
down exact expressions for arbitrary electron operators 
in terms of the bond operators. First, for the electron 
spin operators 
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S2a = 2'^2a'^a6C26, 



(4) 



we have the following expressions which generalize those 
in Ref |^ 

Sic = ^{s^ta + tls - ieap^l^t^) + ^<Tabh\ahlb 



^ Sla + \<bh\ah 
1 

2 



lb 



S2a = + t^S + ieap^t^p^) + ^'yZbhlah2b 



= S2a + -^CTabf^lalT- 



2b, 



(5) 



where eajs-y is the third-rank antisymmetric tensor with 
£xyz = 1- By considering various matrix elements of the 
electron creation operators, we can also obtain 
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ia = h\ad + --^£afcst/l2fc - -^SacC^cbtlh 
4a = hl^d + -^Eafest/lifc + -^SacCTcbt^hib. (6) 
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Indeed, the expression can be obtained (|[) after re- 
peated application of the constraint (^. In a similar 
manner, we can also obtain from (or by direct con- 
sideration of matrix elements) the useful expressions 



'S'lQ5'2ct — '5'iqS'2q 



1 

-1 

4 



c\^Cia = 1 - h\^h2a ~ d) d 
c\aC2a = 1 ^ h\^hia - d'^ d. 



(7) 



We have now completed the definition of the bond op- 
erator formalism. By application of (^,^0) an arbitrary 
Hamiltonian can be written down in terms of the bond 
operators. Notice that the right hand sides of all of these 
equations commute with constraint and so will the 
resulting Hamiltonian. For most systems of interest, the 
Hamiltonian will only contain terms which are quartic 
in the bond operators-we will merely apply the simplest 
possible Hartree-Fock-BCS theory to solve such a model. 

At this point, it is worth pointing an important dis- 
tinction between the present approach and the familiar 
"slave" boson and fermion theories. A crucial advantage 
of our approach is that there are no long-range gauge 
forces associated with the fluctuations about our simple 
Hartree-Fock-BCS theory: the s boson is strongly con- 
densed in all the phases, and so the U(l) gauge sym- 
metry associated with (^) is always badly broken. Con- 
sequently, the quantum numbers of the bare excitations 
are also those of the renormalized quasiparticles; this is 
rarely the case in the slave particle approaches. In par- 
ticular, in the insulating phases, the main elementary 
excitation will be the bosonic ta particle; in the param- 
agnet this is the 5 = 1 exciton, while in the Neel state 
it reduces to the two spin waves. In the superconduct- 
ing states, we will find that the ta excitations persist as 
the S = \ spin exciton, and the /iia, ft.2a quanta turn 
into the fermionic, S — 1/2 Bogoliubov quasiparticles. 
So there is a direct and transparent relationship between 
the quantum numbers of the bond operators and those of 
the elementary excitations of all the ground states. This 
is main benefit of our approach, and we are not aware 
of any previous theory which has satisfied these criteria. 
Of course, some of the quantum numbers of our confining 
theory can also be obtained in a direct Hartree-Fock-BCS 
theory of an extended Hubbard model of the electrons. 
However, this has to be followed by an RPA analysis of 
collective modes to obtain the S = 1 exciton; moreover, 
in such theories, the largest energy scale in the problem, 
the charge gap in the Mott insulator, is unphysically re- 
lated to the strength of the magnetic or charge order, and 
is hence grossly underestimated. 

We conclude this section by noting how the bond op- 
erator theory satisfies criteria (i)-(iii) of Section |. (i) As 
just noted, the ta quanta are gapped S = \ excitons in 
the insulating paramagnet, and persist as sharp excita- 
tions in the superconductor for a finite range of doping. 




(b) 



FIG. 1. Definition of the Hamiltonian J/" in (|8|). The 
two- leg ladder is in (a), while the square lattice is viewed as 
a set of coupled two- leg ladders in (b). The thick links have 



J, v., 



Jij = A J, and Vij — V, while the dashed links have tij = t' , 



Jij = A' J, and Vij — V. 



(ii) Any mean-field theory using the bond operator for- 
malism will prefer the bonds on which the operators re- 
side over the others, and this naturally leads to a broken 
lattice symmetry for symmetric Hamiltonians. Indeed, 
the primary weakness of the bond operator formalism is 
that there is no simple way to restore this symmetry, (iii) 
Placing a Zn/Li impurity means that the partner of one 
site has been removed. The excitations of this site (and 
only this site) therefore cannot be described by the bond 
operators above: instead we need a fermionic spinor, ha, 
to create the S — 1/2 state with one electron, and a spin- 
less boson, b, to represent the hole. At zero doping, the 
ha particle constitutes the free S — 1/2 moment near the 
impurity. This particle will be bound near the impurity 
site for a finite range of doping, and so the moment will 
persist for a while in the superconductor. 



III. RESULTS 

We applied the bond-operator method to the t-J-V 
model defined by the Hamiltonian 



H — 1^ tij (^cjaCja + Cja'^iaj + JijSiaSja 

+ Vrjc\aC,ac],,Cjb " c\aC.,a , (8) 



(io) 



where the sum {ij) extends over nearest neighbor pairs 
on a two-leg ladder or a square lattice, and it is implied 
that all states with two electrons on any site have been 
projected out. The are the electron hopping matrix 
elements, the Jij > are the antiferromagnetic exchange 
interactions, the Vij > are nearest neighbor repulsive 
Coulomb interactions, and /x is the chemical potential. 
The values of the tij , Jy and Vij are indicated in Fig |l| 
for both the two-leg ladder and the square lattice. For 
the two-leg ladder we have tij — t, Jij = J, Vij = V on 
the horizontal links, and tij = t, Jij = AJ, and Vij — V 
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on the vertical links. We view the square lattice as a 
collection of adjacent two-leg ladders: then the couplings 
on each ladder are the same as before, while on the links 
between the ladders we have = t', Jij — A' J, and 
Vij = V. This parameterization is chosen so that in 
both cases the exchange interaction decouples into dis- 
connected pairs of sites ("dimers") at A = A' = 0. Fur- 
ther, at A = A' = 1, 7? has the full symmetry of the 
square lattice. However, as discussed earlier, our mean- 
field theory will continue to have a ground state with the 
symmetry of Fig ^ even at these values of A, A': this im- 
plies the presence of spontaneous bond-centered charge 
order of period 2 in the ground state. 

The calculation proceeds by the substitution of the op- 
erator representations in Section ^ into (||) , followed by 
a Hartree-Fock-BCS treatment of all the quartic terms. 
The procedure is quite lengthy, but the computations are 
quite similar to those in earlier work. Details of the cal- 
culation are presented in Appendices |c| and 

Here we will discuss the results of such a calculation. 
All of the phases can be characterized by a specification of 
the non-zero expectation values of various combinations 
of the bond operators in the ground state. 

All states have the non-zero expectation values 

(s) ^ ; (tit^) ^ 0; (t„t„) ^ 0. (9) 

These non-zero values do not break any physical symme- 
tries of the Hamiltonian, and merely serve to break the 
U(l) gauge symmetry associated with the constraint (||). 
As noted earlier, this is fortunate as no long-range gauge 
forces then appear in the fluctuations about our ground 
states. In addition, all phases will also have non-zero 
expectation values of the operators tJ^^Q, and h\ 2a^i,2a, 
which do not break any symmetries of the Hamiltonian, 
and also commute with the constrains (^). 

At zero and non-zero doping we have find also find 
magnetically ordered states. These are characterized by 
the non-zero expectation value 

(ta) ^ 0; (10) 

so condensation of single bosons (which arc 5* = 1 par- 
ticles) leads to the appearance of magnetic order. This 
also implies that the magnetic ordering transition will 
be described by the field theory a 3-component spin- 
vector order parameter. Further, provided the ordering 
wavevector is not exactly equal to the spacing between 
the points of gapless fermionic excitations, the univer- 
sality class of the transition-,!*, the same as that in the 
0(3) non- linear sigma modclQuJ. This is to be contrasted 
with magnetic ordering transitions in Schwinger boson 
theories, which are associated with the condensation of 
5 = 1/2 particles, possibly interacting with each other 
via gauge forces. 

The singlet superconducting states have the anomalous 
expectation values 

{d)^0 ; eab{hi,2ahi,2b) + 0. (11) 



The spatial pattern of these anomalous condensates de- 
termines the symmetry of the Cooper pair wavefunction, 
and whether there are any nodal fermionic quasiparticles. 

Finally, we also found states with co-existing magnetic 
order and superconductivity, in which both ( p^ ) and ( [Tl| ) 
hold. 

Our main results are summarized in Fig |^. First, let 
us review the results at x = 0, in the insulator. Here the 
calculations are very similar to those already considered 
in Ref |2^. For small A', the ground state is an insulating 
paramagnet: the symmetry of the Hamiltonian and the 
ground state is that of Fig |l]b, and the excitations 
have an energy gap: these excitations constitute a 5 = 1 
collective spin resonance (or a S* = 1 exciton). At a 
critical value of A', the excitation gap vanishes, and long- 
range Neel order sets in. 

Now we turn to non-zero doping. In our Hartree- 
Fock-BCS mean-field theory, we find that superconduc- 
tivity appears at any non-zero doping for small A'. As 
one of us has discussed in Ref. ^ (along with references 
to the earlier literature), this is surely an artifact of 
our approximation at very small x: in the presence of 
long-range Coulomb interactions, Wigner crystal states 
will be present for very small hole concentration. Even 
in our theory without long-range Coulomb interactions, 
long-range charge inhomogeneities can also appear in our 
mean-field frameworkfl; however, we will neglect these 
here for simplicity. As indicated in Fig ^ the supercon- 
ductivity persists in a non-magnetic state even at A' = 1 
for X > 0.082. In our present calculation, this small 
X superconducting state has bond-centered charge order 
of period 2, but more complejj-bond-centered charge- 
ordered states are also possiblcDa. We also expect that 
the superconductivity will co-exist with antiferromag- 
netism, but we did not undertake a complete solution 
of the mean-field equations within the antiferromagnetic 
phase because of the complexity of the calculation. As in 
the insulator, the onset of antiferromagnetim was deter- 
mined by the point where the ta boson had a vanishing 
gap in its spectrum. 

We characterized the superconducting state by deter- 
mining pairing amplitudes, A^, j,, which characterize the 
pairing of the holes, /ii,2j in the x and y directions. These 
quantities are described more precisely in Appendices ^ 
and 1^. Their values are plotted as a function of the hole 
concentration, x, in Fig ^j. Note that A^, and Ay are 
roughly of the same magnitude: this means that the su- 
perconductivity is genuinely two-dimensional, and there 
is no regime in which a quasi-onc-dimensional Luther- 
Emery liquid-like behavior holds. Indeed, it is even pos- 
sible for A^; to be larger than Ay, which would make 
the superconductivity stronger in the x direction rather 
than in the vertical "stripe" direction. This latter phe- 
nomenon appears to be similar in spirit to the ananlyses 
of Refs |6|-||. 

For small x, the pairing amplitudes are non-zero but 
quite small; moreover, the pairing amplitudes in the hor- 
izontal and vertical directions have the same sign in our 
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FIG. 2. Phase diagram of the square lattice Hamiltonian 
H in (^) for the configuration shown in Fig Results are 
shown here for t' = t, \ = 1 and t/J = 3.0 as a function 
of the hole concentration x and 1 — A'. The model has full 
square lattice symmetry at 1 — A' = 0, and the exchange in- 
teraction separates into that on decoupled two-leg ladders at 
1 — A' = 1. The nearest-neighbor Coulomb repulsion, V is 
taken to have a very large positive value, and the results are 
indistinguishable from those a.tV = oo. The phase boundary 
towards antiferromagnetism was determined by the point at 
the ta excitation spectrum had a vanishing spin gap, signaling 
the onset of a phase in which (|lo|) was valid. The "weak" and 
"strong" superconductivity distinction is qualitative, and is 
indicated by the rapid change in the physical superconduct- 
ing order parameter in Fig H around x = 0.2. The antifer- 
romagnetic phase co-exists with weak superconductivity at 
non-zero x. The nearly vertical phase boundary indicates the 
quantum phase transition at which nodal fermions first ap- 
pear: this transition is discussed in Section IV C denotes 
the point where the spontaneous bond-centered charge order 
of the square lattice is expected to disappear with increas- 
ing x; the bond-centered charge order is explicitly present in 
the Hamiltonian for 1 — A' > 0, and is spontaneous only for 
non-zero x before the point C and 1 — A' = 0. The present 
bond operator approach was not used to determine the posi- 
tion of C; nevertheless, following Ref ^ we are able topresent 
a theory of the critical properties of C in Appendix 
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FIG. 3. Pairing amplitudes (A) as a function of hole con- 
centration (x) when t/J — 3.0, and A = A — 1.0, i.e. the 
isotropic case. The pairing amplitude in the vertical direc- 
tion of Fig 1^ is Ay, and that in the horizontal direction is 
A^. For X < 0.082, there is instability to long-range antifer- 
romagnetism due to the collapse of the gap of S = 1 particle 
excitation, ta- When 0.082 < a; < 0.192, non- magnetic su- 
perconducting state is stable, and the pairing amplitudes are 
of the same sign, but of very small value. For doping concen- 
tration larger than roughly 0.2, the pairing amplitudes have 
opposite signs, which, in turn, gives rise to nodal fermions. 



mean field theory. We will refer to this as the weak super- 
conductivity regime. It is possible that, upon including 
quantum fluctuations, the ground state in this weak pair- 
ing regime is easily susceptible to a quantum transition 
to an insulating state with some additional translational 
symmetry breaking. Also note that the pairing amplitude 
is nonzero at the boundary of the onset of antiferromag- 
netic order: we therefore expect the superconductivity to 
survive within the antiferromagnetic phase. 

For larger x, the pairing amplitude in Fig ^ in- 
creases rapidly and we reach a "strong superconductiv- 
ity" regime. The pairing amplitudes now have opposite 
signs in the horizontal and vertical directions, and this 
permits gapless nodal fermionic quasiparticle excitations, 
as shown explicitly in our results below. However, there 
is one important feature associated with the appearance 
of the nodal particles that is worth emphasizing. Origi- 
nally, our calculations were carried out with the nearest- 
neighbor Coulomb repulsion V = 0. In this case we found 
that the pairing of the holes occurred primarily by the 
condensation of the d bosons. Because of the resulting 
very short-range pairing we found no nodal points, even 
though the pairing amplitude had opposite signs in the 
horizontal and vertical directions. It was only when we 
had turned on a V of the order of the bandwidth, which 
significantly reduced the amplitude of the d boson con- 
densate and made the hole pairing more long-ranged, did 
we find the appearance of the nodal excitations. 
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FIG. 4. Dispersion curves of the bonding and anti-bonding 
= 1/2 fermionic excitations in the weak pairing regime with 
X — 0.1, t/J = 3.0, and A = A = t /t = 1. There appear 
to be fermionic excitations at energies very close to 0, but 
the blow-up in Fig ^ shows that there is indeed a gap in the 
fermionic spectrum. The x-a:xis is chosen to take a represen- 
tative straight path connecting points in the first Brillouin 
zone: from {px,Py) = (0,0) to (0, vr) to {n/2,n) to {n/2,0) to 
(0,0). 



We now describe our results for the elementary excita- 
tions of the superconducting states in our phase diagram. 
First, we discuss the fermionic 8=1/2 excitations; in 
principle, these are observable in photoemission or tun- 
neling experiments. Note that the unit cell of Fig ^ 
has two sites, and so the first Brillouin zone extends be- 
tween — 7r/2 and 7r/2 in the x direction, and between — tt 
and TT (as usual) in the y direction. Consequently, half 
the fermionic excitations in the first Brillouin zone of the 
original square lattice will be folded back, by a Bragg re- 
flection in the vertical planes at ±7r/2, into a second band 
in the first Brillouin zone of the lattice of Fig |l|b. We 
will refer to these two bands as the 'bonding' and 'anti- 
bonding' bands, based upon their wavefunctions within 
the dimers, as discussed in Appendix ^ We show the 
S = 1/2 fermionic bands in the weak-pairing regime in 
Figs I and |. There is a gap across the entire Brillouin 
zone. However, because of the very small pairing ampli- 
tude, the fermionic excitation energy becomes very small 
along an incipient "Fermi surface" in the Brillouin zone. 
This Fermi surface line is indicated by the solid line in the 
inset of Fig ||. The position of this Fermi surface is sim- 
ilar to that in the computation by SushkovO, who also 
discussed its relationship to photoemission experiments. 

Next, we turn to a discussion of the fermionic excita- 
tions at larger x, where the superconductivity is stronger. 
The 5 = 1/2 excitation spectra are now shown in Figs ^ 
and 1^. Overall, the results at higher energies, are quite 
similar to those in Figs ^ and ^, but there is now a dra- 
matic difference at lower energies. There is a gapless 




(0,0) 



(71/2,71) 
(Px'Py) 

FIG. 5. Blow-up of the fermionic dispersion in Fig ^ at the 
lowest non-zero energies. Notice that there is indeed a finite 
gap in the fermionic spectrum. This is also made clear in the 
inset: the solid line shows the line where the fermion kinetic 
energy vanishes (the incipient Fermi surface) while the dashed 
line indicates where the fermion pairing amplitude vanishes. 
For gapless excitations, both quantities have to vanish, and 
the absence of a crossing point between the lines indicates 
that a gap is always present. The solid line in the inset is also 
the location of the incipient Fermi surface which has been 
quenched by pairing — the position of this surface is similar to 
that in Ref |l[ 
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FIG. 6. Dispersion curves of the bonding and anti-bonding 
5 = 1/2 fermionic excitations with x = 0.3, t/J — 3.0, and 
A = A — t /t = 1 (as in Fig ^ but at a larger doping). At 
this value of x, there is a nodal quasiparticle excitation not 
too far from (n /2,-k /2). 
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FIG. 7. Blow-up of the fermionic dispersion in Fig ^ at the 
lowest non-zero energies (as in Fig ^ but at a larger doping) . 
Notice that there is clearly a gapless nodal point not too far 
from (7r/2, 7r/2), as is also clear from the inset, where the two 
curves now cross. The pairing gap near the "antinodal point" , 
(tt, 0), remains finite. 



nodal-point near (7r/2, 7r/2), and this is made possible by 
the d- wave-like pairing amplitudes. The spectra in Fig 
and 1^ cannot be continuously connected, and there must 
be a singular phase transition point at which the nodal 
fermion first appears. The position of this phase bound- 
ary is indicated in Fig ^ and the nature of the phase 
transition will be discussed below in Section 

Finally, wc consider the S — 1 excitations, which are 
visible in neutron scattering experiments, There are two 
categories of such excitations. First, we have the S = 1, 
ta particles, with a definite energy-momentum relation: 
these may be viewed as a collective spin mode (or a triplet 
exciton) which goes soft at the antiferromagnetic order- 
ing transition. These are the only S — 1 excitations 
in the undoped antiferromagnet, but, connect smoothly 
to corresponding excitations in the doped antiferromag- 
net. The second class of S" = 1 excitations are scattering 
states of the fermionic S = 1/2 excitations just discussed. 
These do not have a definite energy momentum relation, 
but exist over a continuum in a range of energies at any 
given momentum. We show a plot of the spectra of these 
two S — 1 excitations in Figs || and ^. 

At low doping, in Fig |8[ not that the ta excitation 
has a well- formed minimum at (0,7r), and this minimum 
is below the lower-bound of the two-particle continuum. 
So the ta excitations are absolutely stable towards de- 
cay into 2 fermionic 5=1/2 excitations over this range 
of energy and momenta. This stable 5 = 1 excitonic 
excitation will lead to a resonance peak (in principle, in- 
finitely sharp) in the neutron scattering cross-section at 
the energy and momentum of the ta particle. Also, as 
the doping is lowered, the minimum of the ta disper- 
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lower bound of two-particle excitation continuum 




^'%,0) (0,7t) (71/2,71) (71/2,0) (0,0) 

(Px'Py) 

FIG. 8. Dispersion curve (full line) of 5 = 1, ta particle 
excitation with x = 0.1, t/J = 3.0, and A = A = t /t = 1. 
The pairing is weak for this value of doping. Also shown 
(dashed line) is the lower bound of the two-particle continuum 
made up of a pair of 5* = 1/2 excitations from Fig ^ and ^. 
Note that the minimum of the ta excitation is at (0,7r), but 
because of the halving of the Brillouin zone this wavevector 
is crystalographically equivalent to (7r,7r). 
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FIG. 9. As in Fig ^, but in the strong pairing region with 
nodal quasiparticles at a; = 0.3. 
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sion gets lowered, until it hits zero also in a region which 
has no overlap with the fermionic two-particle contin- 
uum. This signals the transition to the antiferromagnetic 
phase: the absence of low energy fermionic excitations 
at near the minimum ta energy implies that the criti- 
cal theory for the transition remains identical to that for 
the corresponding transition iPirthc insulator, which has 
been discussed in much detailocj. The condensation of 
the ta bosons occurs at momentum (0, tt), but note that 
because of the doubling of the unit cell this momentum 
is equivalent to (7r,7r). Determination of the microscopic 
configuration of the spins using the operator representa- 
tion in Section || shows that the mean moments indeed 
have the staggered arrangement associated with the con- 
ventional, two-sublattice Neel state. 

The higher doping S = 1 excitations are sketched in 
Fig ||. Now the ta excitations are well within the two- 
fermionic-particle continuum, and so will quickly decay 
and lose their identity. Note that the doping at which 
the ta excitations ceases to be stable at any momentum 
does not (in general) co-incide with the point at which 
nodal fermionic quasiparticles appear. As long as the 
nodal points are not exactly at (7r/2,7r/2), the ta exci- 
tations can be stable at very low energies near (7r,7r). 
In the present computation, the nodal points are quite 
close to (7r/2,7r/2) and to the ta excitation is quickly 
quenched. In a more general model, with second neigh- 
bor hopping, the ta excitations could be stable even in 
an isotropic d-wave superconductor, provided we moved 
the nodal points sufficiently far from {tt/2, tt/2). 



IV. ONSET OF GAPLESS NODAL FERMIONIC 
EXCITATIONS 

The phase diagram of Fig || shows a phase bound- 
ary demarking superconductors with and without nodal 
quasiparticles. This is a quantum phase transition in the 
sense that there is a (weak) non-analyticity in the ground 
state energy as a function of doping at this point. We will 
present a simplified treatment which captures its essen- 
tial universal features. Such a transitkm was discussed 
earlier in Ref and by Granath et aZ.c£l; the latter au- 
thors also reached conclusions on the universal properties 
which agree with our discussion here. A related, but dis- 
tinct, theory for the annihilation of nodal particles in 
a d-wave supercpaductor was also discussed by Duncan 
and Sa de MelocZl. They considered an isotropic super- 
conductor as a function of electron density, and found 
that the nodal points vanished when all four of them col- 
lided ai k — 0. This transition is differs from that in 
the anisotropic superconductors of interest here, where 
the nodal points collide only in pairs, and the quantum 
critical points belong to distinct universality classes. 

We approach this transition from the side of the d- 
wave superconductor with full square lattice symmetry. 
Here, the fermionic excitations are described by the BCS 
Hamiltonain 



H 



BCS 



E 

k 



A(fc) 



-kb 



H. 



(12) 



where, in the simplest nearest-neighbor model — 
—2t{coskx + cos ky) — fi, and the pairing energy A(A;) = 
Ao(cos fc^— cos ky). We now assume that there is an onset 
of bond-centered charge order at wavevector G = (tt/o, 0) 
of amplitude ■0sp — because this charge order wave is in 
the X direction, ipsp is the real part of the more general 
order parameter, considered in Section | and Ap- 
pendix This order will lead to a modulation in the 
fermion hopping matrix element and the pairing inter- 
action at the wavevector G, and so induce the following 
additional terms in the Hamiltonian for the fermionic ex- 
citations 



Hsp = itpsp ^ a{k)cl^^Ck+G,a 



Kk) 

2 



C'-k+G.b + H. c. 



, (13) 



where a(k) = wisinkx and b{k) — W2siiika:, with ^11^2 
some constants. These last factors of sin kx are a conper 
quence of the bond-centered nature of the charge ordercS, 
but the results of this section are not crucially depen- 
dent upon this fact; similar results will apply also to 
site-centered charge orders. 

It is possible to diagonalize the Hamitonian Hbcs + 
Hsp and determine the fermionic excitation spectrum of 
the state with co-existing superconductivity and charge 
order. The energy eigenvalues are 



Ai ± (A? 



^^2 



Al 



1/2 



1/2 



(14) 



where 



Ai EE [e'ik) -f e^k + G) + A2(fc) + A''{k -f G)] /2 

+ ^%[a'{k)+h\k)\ 
A2 = A(fc)e(fc + G) + A(fc + G)e{k) - 2tP^pa{k)b{k) 
A3^<[a2(fc)-62(A:)] 

+ A{k)A{k + G) - e{k)e{k + G). 

It is instructive to examine the evolution of the zeros 
of ([ij) as a function of ipsp- For V^gp = 0, there are four 
symmetric nodal points determined by the solutions of 
e(fc) = and A(fc) = 0. It is assumed (as is the case at 
non-zero doping and in the absence of particle-hole sym- 
metry) , that the wavevector separating any two of these 
nodal points is not equal to G. If it was equal to G, 
then the nodal points would be gapped at an infinitesi- 
mal value of ^ps■p^ and a non-trivial theory would describe 
the quantum critical fluctuations, as has already been 
discussed in Ref ^ Turning to the more general situa- 
tion, where the separation between the nodal points is 
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not equal to G, the quantum critical theory describing 
the onset of a non-zero 4'sp is reviewed in Appendix ^ — 
the leading critical singularities do not involve the nodal 
fermions. As V'sp increases, examination of (^) shows 
that the four nodal points move towards the new Bril- 
louin zone boundary at — ±7r/2. Eventually, these 
points collide in pairs at 



Po = (7r/2,poa;) 



(15) 



(and symmetry related points) and ipsp = V'spi there 
are no nodal points for ■0sp > V'sp- We are interested here 
in describing the quantum critical theory of this nodal 
point collision. Clearly, the fluctuations of ipsp will not 
be critical at this point, as the transition occurs at a 
non-zero value of '0sp. 

First, let us determine the values of ip^^ and po. Using 
the facts that e(po) = e(Po + G) and A(po) = A(po + G), 
the condition for the presence of gapless nodal points in 
(nj) becomes 




e(Po) = ■>Pspa{po 



A(po) = ^spHpo) 



(16) 



The solution of these two equations determines the two 
unknowns Vsp ^^'^ Pox- At this same point, the sec- 
ond eigenvalue in (|lj) remains non-zero and finite. This 
fermionic mode will play no role in the critical theory, and 
so it pays to perform a canonical transformation at an 
early stage to project it out. To leading order in the dis- 
tance from the critical point, this merely means that we 
have to take only the linear combination oicka and Ck+c.a 
which appears in gapless eigenvalue associated with (16). 



So, we introduce a new fermionic degree of freedom, 7|^, 
where the small momentum q is measured as a deviation 
from Po (and a corresponding fermionic mode which 
resides at momenta near — po); using the structure of the 
gapless eigenvalue at V'sp — V'sp ^^"^ = po, we see that 
we should parameterize 



Ck+G,a = i£abf, 



qb 



(17) 



where k = pq + q, and q is small; a similar parameteri- 
zation is made near —po with /~ . Finally, we insert (|l7| ) 
into Hbcs + ^spj and expand in V-'sp ~ V-'sp ^^.'^ i'^ S^a- 
dients of This leads to the following effective action 
for the critical theory 



5*/ = / (frdr 



d 



d 



1 92 



Sabfa \ ^2-Q^ 



dr 

1 



dy 2mi dx^ 



92 



2m2 dx^ 



+ 52 /+ + H.c. 



(18) 



where r is imaginary time, r = (x, y) is the spatial co- 
ordinate, 5i,2 = Wl,2(V'sp - V'sp); and Ui,2, "11,2, w[ 2 

are constants dependent upon the detailed momentum 
dependence of e{k) and A(A:). The fermionic cigencnergy 
of (na) is easily determined; it is 



FIG. 10. Evolution of the gapless nodal points in the 
fermionic excitation spectrum as a function of V'sp- The 
square contains the first Brillouin zone of the original square 
lattice, extending between kx^y — ±7r. The dark circles are 
the positions of the nodal points in this d-wave supercon- 
ducvtor. Onset of a non-zero Vsp introduces Bragg reflection 
planes at kx = ±7r/2 and images of the nodal points indicated 
by the grey circles. As V'sp increases, the nodal points move 
towards the Bragg-|refiection planes and annihilate each other 
when they coUideu. The critical theory of this transition is 
discussed in Section 



5i + viqy 



2m 1 



V2qy 



2m2 



1/2 



(19) 



The positions of the nodal points, if present, are easily 



determined from (|19|). Assuming wi,2, 'ml^2^ 



are all 



positive, then for w'l/vi > w'2/v2 and mivi < m2V2, 
nodal points are present for V'sp < V'sp, but not for V'sp > 



V'sp (similar results hold for other signs and magnitudes 



of the various coupling constants). The trajectory of the 
nodal points is sketched schematically in Fig 10: they 
move on a parabolic trajectory with \qy\ ~ before 
colliding along the Bragg reflection planes at kx — ±tt/2. 

The availability of the action Sf also allows one to 
determine the consequences of the interactions near the 
quantum critical point. Notice that at the critical point, 
^1 =1^2 = 0, Sf is invariant under the scale transforma- 
tion T —> t/s, y y /s, X x/^/s, and / s^^'^f- The 
simplest allowed four fermion coupling is ~ (/V)^, and 
power-counting shows that its co-efficient has a negative 
scaling dimension of —1/2. So this interaction is an irrel- 
evant perturbation at the quantum critical pointCJ. In a 
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similar manner, it is not difficult to show that all other 
perturbations of Sf are irrelevant, and so Sf is the com- 
plete critical theory of the transition involving the onset 
of nodal fermionic excitations. 



V. COUPLING OF SPIN PEIERLS ORDER TO 
PHONONS 

A central actor in all the considerations so far has 
been the spin-Peierls order parameter, ^'sp. Experimen- 
tal observations of the fluctuations of this order parame- 
ter would certainly be helpful in resolving the theoretical 
issues. However, this is a spin-singlet, charge zero mode 
and is only observable through its couplings to the ionic 
displacements. This section will therefore present a sim- 
plified discussion of the coupling between the dynamics of 
the spin Peierls order parameter and the phonon modes of 
the Cu02 plane. Indeed, it is probable that the phonons 
are more than merely spectators of the spin dynamics, 
and the spin-phonon coupling may play an important 
role in selecting between different charge orderings and 
in influencing the nature of the electronic ground state: 
the important phonon and spin excitation energies are 
roughly of the same order, and so a coupled dynamical 
model will be necessary for a complete microscopic un- 
derstanding. 

The spin-Peierls order parameter couples most directly 
to the "Peierls-active" phonon mode sketched as mode C 
m Fig |ll|: this is a staggered displacement of the Cu 
ions at the wavevector (tt, 0). Notice that the O ions 
are stationary in this phonon mode. As we will describe 
below, the coupling between this Peierls-active phonon 
and the spin-Peierls order parameter is likely to lead to 
an additional low energy peak in the dynamic phonon 
structure factor at the energy scale of the characteris- 
tic spin-Pcicrls fluctuations. However, neutron-scattering 
studies of the phonon spectra have so faiOcjEil focused 
on optical phonons involving motion primarily of the O 
ions (the mode between A and B in Fig [ill) : the smaller 
O ion mass, Mq ~ Mcu/4, increases the displacement 
of the O ions and makes these easier to observe. The 
O ion phonon at wavevector (tt, 0) (mode B in Fig 
is not Peierls-active in the sense just noted. |-J3owever, 
following the suggestion of McQueeney et aZ.Eil, we will 
describe below a non-linear coupling in the spin-phonon 
Hamiltonian which leads to signatures of the spin-Pcicrls 
fluctuations in the A-B phonon mode at the wavevector 
(V2,0). 

It is also worth—mentioning the recent work of Khal- 
iuUin and HorschEj. In contrast to our focus on bond- 
centered charge fluctuations, they considered Cu site- 
centered charge density fluctuations at wavevectors near 
(tt, 0) (in the notation of NayakH, we are examining px 
density wave correlations at (tt, 0), while Ref |5^ discusses 
s density wave correlations at (tt, 0)). These Cu site- 
centered fluctuations do indeed couple linearly to the 




FIG. 11. Phonon frequencies (vertical axis) of the 
one-dimensional Cu-0 chain. The frequency of the mode A 
is [2K{Mcu + Mo)/(McuMo)]'/^ that of B is [2A'/Mo]''", 
and that of C is [2K/McuY^'^. Only mode C is "Peierls ac- 
tive" i.e. couples linearly to the order parameter 'I'sp, which 
measures the amplitude of a bond-centered charge order wave 
at wavevector (7r/a,0). 
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phonon mode in the vicinity of the point B in Fig |T|: 
this led to a broadening and softening of the phonon near 
B. As we will see below, our bond-centered spin-Peierls 
correlations instead modify the A-B phonon mode near 
(7r/2, 0) in a manner which is consistent with experimen- 
tal observations. 

We will work here with a simplified one-dimensional 
model of the phonons in the Cu02 plane. All the phonons 
we are interested in are polarized both in wavevector and 
ionic displacement along the (1, 0) axis, and so neglecting 
the second dimension is not too serious. We will consider 
only a single chain of alternating Cu and O ions, and con- 
sider displacements of the ions along the chain direction 
(see Fig 11). A more complete treatment is certainly 
possible, but we will not attempt it here: we hope that 
with improved experimental resolution, a more precise 
and microscopic theoretical study of the two-dimensional 
phonon modes will be carried out, along the lines of the 
analyses of Refs 53 for CuGeOa. 

We begin by specifying our toy one-dimensional model 
of the phonons in an alternating chain of Cu and O ions. 
We place the Cu ions on the sites, i (integer), of a chain 
with spacing a, while the O ions are at the centers of the 
links of the chain. We denote the displacement of these 
ions by m and Vi respectively, where the i'th O ion is 
taken to be to immediately to the right of the i'th Cu 
ion. The harmonic action for the phonon modes is (in 
imaginary time, r): 



5'ph ^ J drY^ 



Mcu ( du 



dr 



Mq / dwj 
~2~ \ d7 



(20) 



where X is a "spring constant" which determines the 
phonon frequencies. It is a simple matter to diagonal- 
ize 5ph and obtain the phonon normal mode frequencies: 
they are 



K 



McuMc 



Mcu + Mo ± ((Mcu - Mc 



-4McuA^ocos2(fc/2)) 



1/2 



(21) 



and are sketched in Fig |T|. 

We also have to consider the dynamics of the spin 
Peierls order introduced in Section | and consid- 
ered in Section |^ and Appendix ^ In our present 
one-dimensional toy model, we need only consider Tpsp = 
Re[5'sp]. For this order parameter ipsv ^ th e dynamics of 
an interacting effective action like S* in ( Al ) is assumed to 
be captured by the following effective quadratic actionQ: 



T 



(22) 



X, 



sp 



+ 2a^c\{l - cos(fc)) -f -f rsp|cj„|. 



Here ujn is a Matsubara frequency, and Agp and Fgp are 
effective energy scales determining the mean frequency 



and damping of the '0sp fluctuations; the values of Agp 
an d Fg p are determined by the non-linear interactions 
inJAl). Also we have replaced the spatial gradient in 
(|A1) by a nearest-neighbor lattice derivative in (^2|). As 
one approaches the onset of spin-Peierls order (approach- 
ing point C in Fig |^ from the right), the value of Agp 
will decrease to zero, while Fgp becomes of order T, rep- 
resenting the damping of the order parameter mode in 
the quantum-critical region. We will work here in the 
Asp > regime, staying to the right of C in Fig ^. The 
velocity ci should be order the spin-wave velocity, and 
this is about 10 times larger than the velocity of the 
acoustic phonon mode, C, in Fig |ll]. 

Finally, wc_h.a|jie to couple il^s-p to the phonon modes. 
This couplingp3L3 arises from the dependence of the ex- 
change constant J between neighbor Cu spins on the 
displacements of the Cu and O ions. We assume that 
J ^ tpj^, where tpd, the overlap between neighboring O 
and Cu orbitals is a function of {ui — Vi) and {ui — Vi-i). 
Expanding J in derivatives of these variables, we obtain 
first the simple linear coupling 



Sic 



-l)Vsp,i (wi+1 



0] 



(23) 



where A is the linear coupling constant, and V'sp.i is the 
spin-Peierls order parameter in real space; this naturally 
resides on the centers of the bonds, and we locate V'sp,i on 
the same O site as Vi . As we noted earlier, Sic couples 
the spin-Peierls order most strongly to the phonons in 
the vicinity of the point C in Fig |ll|, and will lead to its 
broadening and softening. We are interested here primar- 
ily in the modifications of the optical phonon mode A-B, 
and so we will neglect Sic in our computations below. 

As indicated at the beginning of this section, the im- 
portant effect on the A-B optical phonon arises from 
a non-linear spin-phonon coupling . Expanding to one 
higher order in the phonon displacements, the coupling 
between ijjs-p and the phonon displacements can be writ- 
ten in the form 



S2c = J dr'Y^ [(-l)Vsp,i jli ("i+i - u^f 



72 (u 



i+1 



2v. 



>-} 



(24) 



where 71.2 are the non-linear spin phonon coupling con- 
stants. 

We examined the properties of Sph + Ssp + S2C in a sim- 
ple one-loop approximation: we neglect the lu^ phonon in 
( pT| ) , computed the self-energy of the w+ optical phonon 
at order 7^ 2 • Finally, to compare to neutron scattering 
experiments, we computed 



D{k,u) = { bcuuik,uj) + bovik,Lo)e'''^^ 



(25) 



where 5cu = 7.718 and bo — 5.803 are the, neutron scat- 
tering lengths of Cu and O respectively!^. The results 



12 





Asp/wo = 0.15 




Asp/o)o = 0.10 




FIG. 12. Plot of the phonon spectral density lmD{k, uj) de- 
fined in (p5|). The figure shows the influence of the spin-Peierls 
fluctuations described by (^2|) on the optical phonon mode 
A-B in Fig the two degrees of freedom are coupled 
by the non- linear terms in (Bll). We used the parameters 
LOO = ^if/(l/Mcu + 1/Mo) r,p = 0.5a;o, T = O.Iljo, 
c = 10.0a;o, and 71 = 72 = CSifwo- The figures show the 
evolution in the spectrum as a function of Asp/a;o. 
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FIG. 13. As in Fig y but for Fsp = 0.2ujo, T = O.lojo, 
c = S.OciJo, and 71 = 72 = Q.^Kuq. 
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are shown in Figs ^ and for a representative set 
of parameters. As Agp decreases, there is initially a 
broadening and a sharpening of the phonon mode near 
k = 7r/2. For even smaller Agp we see the incipient 
Brillouin zone boundary of the doubled unit cell devel- 
oping near k = tt/2: this is reflected in the apparent 
discontinuity of the dispersion and the appearance of 
"shadow" phonon bands reflected across the incipient 
Brillouin zone boundary. The first of these features is 
in excellent accocd-with all the available neutron scatter- 
ing experimentsEilcj'EJ, and there has even been a claim 
of the observation of shadow bandsu, although this has 
not been confirmed. 

Another interesting feature of Figs |l^ and |l^ is that 
the intensity of the phonon is larger near the Brillouin 
zone boundary. This is a 'form factor' effect, and arises 
primarily from the interplay of the different scattering 
lengths in (|2^), and the distribution of the phonon nor- 
mal mode between the O and Cu sites. Curiously, pre- 
cisely such a dominance pf.pLDtensity near fc = tt is ob- 
served in the experimentsoEiJ. 



VI. CONCLUSIONS 

This paper has introduced a bond operator formalism 
to represent the degrees of freedom of doped antiferro- 
magnets. While the approach bears some similarities to 
the popular "slave" boson and fermion techniques, the 
implementation leads to theories with a rather different 
structure: there are no long-range gauge forces in the 
fluctuations about any reasonable saddle-point in the lat- 
ter formalism, and the quantum numbers of the true ele- 
mentary excitations are simply connected to those of the 
microscopic operators. This is a powerful advantage of 
the bond operator method, and allows much information 
to be gleaned from simple mean-field Hatree-Fock-BCS 
computations. The main disadvantage of the method is 
that it requires a pairing of the sites into bonds at the 
outset. Such a pairing is naturally present in systems 
with spontaneous or explicit bond-centered charge order 
(as in a spin-Peierls states), and it is for these systems 
that the approach is best suited. 

The defining equations of the bond operator formal- 
ism were presented in Section ||, and are contained in 
Eqns (|l[]^). Next, we applied this formalism to the two- 
dimensional doped antiferromagnet sketched in Fig nb. 
The main phase diagram of the model is sketched in Fig 2, 
and the important properties of the phases are summa- 
rized in its caption. The remaining figures in Fig |l| and 
their captions summarizing the excitation spectra of the 
phases — the fermionic, S = 1/2, spectra are in Figs H, 
^, 1^, and 1^, while the 5 = 1 excitations are in Figs^ 
and H; the latter consist of a bosonic, S = 1 exciton, 
and the two-particle continua of the fermionic 5 = 1/2 
excitations. 

The connection of these result to recent neutron scat- 



tering measurements of phonon spectraEJ 
sidered in Section 0. This section did not use the bond 
operator formalism. Instead, it considered the conse- 
quences of incipient bond-centered charge order in an 
isotropic state, in a simple pcrturbative calculation using 
a spin-phonon model. We introduced an order parame- 
ter, ^'sp, characterizing the ordering pattern associated 
with FiB-l|b, and wrote down a phenomenological free 
energyLrE3describing its fluctuations in spacetime in the 
vicinity of the quantum critical point labeled by the point 
C in Fig 0. On the symmetric side of C (the region with 
(^sp) = 0), these fluctuations were controlled by an en- 
ergy scale, Asp, which vanished as C was approached, 
and we considered the evolution of the optical phonon 
spectra as a function of decreasing Agpi these results are 
contained in Figs |l^ and |l^. j— . 

In the unifying language introduced by NayakCj, the 
bond-centered charge order parameter ^'gp may be con- 
sidered as the amplitude of a Px density wave at wavevec- 
tor (tt, 0). The ordinary s density wave at wavevector 
(tt, 0) is associated with Cu site-centered charge order, 
and these may also be strong in the lightly doped an- 
tiferromagnet, especially in the region with long-range 
magnetic order. Their influence on the optical phonon 
spectra was considered in Ref ^2| — the primary effect was 
a broadening of the O optical phonon near (7r,0). This 
should be contrasted with the influence of the 'I'sp fluc- 
tuations described above — the strongest effect was near 
(7r/2, 0) where the phonon dispersion sharpened consider- 
ably, along with a significant amount of broadening. The 
latter effects are_.clearly seen in recent neutron scatter- 
ing experimeatsEjO, although some softening at (tt, 0) is 
also apparenttil. We hope that higher precision and more 
detailed neutron scattering experiments will be under- 
taken, and along with more microscopic theoretical com- 
putations, these should help sort out the relative roles 
of site- and bond- centered charge order as a function of 
increasing doping. 

In addition to fluctuating charge-order modes detected 
in phonon scattering, it would also be useful to study sys- 
tems in which the charge order is required to be static; 
in such situations, atomic resolution STM studies should 
yield much useful information on the microstructure of 
the charge order. Our physical picture implies that static 
charge order should be present in situations in which 
both magnetic and superconducting order have been sup- 
pressed (systems with one of these orders may only have 
fluctuating charge order). A convenient way to achieve 
this is by application of a strong magnetic field on un- 
derdoped samplesE3. A phenomenological theory of the 
phase diagram in a magnetic field has been provided re- 
cently in Ref the "normal" state in this phase di- 
agram is a very attractive candidate to bond-centered 
charge order. It would be especially interesting to con- 
duct STM measurements on the strongly underdopied 
YBCO crystals that have become available recentl5iO, 
after superconductivity has been suppressed by a static 
magnetic field. An alternative is to look for charge or- 
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der in STM studies in which the superconductivity has 
only been locally suppressed, as is the cas^n the cores 
of vortices in the superconducting ordertllEHI. However, 
the short-range nature of the suppression means that 
charge order is not required to appear, and may remain 
dynamic — this piakes this approach less attractive. Re- 
cent indicationaSj of mesoscale self-segregation of charge 
carriers in bulk samples also naturally raise the possi- 
bility of bond charge order in the lower density regions 
which (presumably) have suppressed superconductivity. 



A. Review of related work 

Before closing this paper, we present brief critical re- 
views of some related work in the literature, and relate 
it our results. 

Sushkov and collaborators have studied aspects of 
Mott insulators both at zero and finite doping. At zero 
doping, they have examined the square lattice 5=1/2 
Heisenbserg antiferromagnet with first (Ji) and second 
(J2) neighbor exchangeEil: with increasing J2/J1 they 
find a transition from the familiar Neel state to a para- 
magnetic state with bond-centered charge order as in 
Fig albeit with a small window of parameters at 
which the two orders co-exist; such a co-existence region 
is allowed in the framework of the earlier field-theoretic 
analysisO. At somewhat larger J2/J1 they find further 
charge ordering which breaks an additional Z2 lattice 
symmetry: this pattern of charge orderijag has also been 
seen in recent DMRG studiesE3. Sushkovcll has also stud- 
ied the evolution of the ground state as a function of 
doping: by a series expansion method he finds that after 
destruction of the Neel order, the ground state has the 
columnar spin-Peierls order of Fig Only at a signifi- 
cantly larger coupling is the full square lattice symmetry 
restored. All of these results are fully consistent with the 
approach outlined in our paper. _ _ 

Mazumdar, Clay, and CampbellE3~E3 have stud- 
ied two-dimensional correlated electron models asso- 
ciated with the organic charge-transfer solids tetra- 
methyl-tetrathiafulvalene (TMTTF), tetramethyl-tetra- 
selenafulvalene (TMTSF), bisethylencdithio-tetrathiaful- 
valene (BEDT-TTF) and bisethylenedithio-tetraselena- 
fulvalene (BETS). They have argued that these materials 
are characterized by quarter-filled band, and produced 
evidence for "bond order wave" in their models at quar- 
ter filling. These states are closely related to the period 4 
bond-centered charge density wave state studied in Ref |^, 
and so also to the insulating and superconducting states 
with the symmetry of Fig Hb stujiied here. 

Kivelson, Fradkin and Emeryl£3 have studied symme- 
try breaking in two-dimensional electronic systems us- 
ing a rather general classification in terms of electronic 
liquid crystal phases. Their formulation (and their ear- 
lier theoryliy of the "spin gap proximity effect") does 
not distinguish between the nature of bond- and site- 



centered charge order that we are paying particular at- 
tention to here. A significant portion of the supercon- 
ducting state in their phase diagram overlaps with a re- 
gion of nematic order. We believe that the orthorhombic 
anisotropy of this nematic order is ultimately driven by 
singlet bond correlations associated with the spin-Peierls 
order in Fig |l|b; consequently this nematic region should 
show bond-centered smectic correlations with a period 
of 2 lattice spacings, as has been claimed in the experi- 
ments of Ref 21, We make this distinction more precise 
by discussing how the spin-Peierls and nematic order|S-are 
coupled. We can characterize the spin-Peierls ordeiO in 
Fig |l|b by a complex Z4 order parameter which takes 
the values 1, i, — 1, — z on the four states obtained by ro- 
tating Fig Pd about any lattice site successively by 90°. 
Similarly, we can introduce a real Ising nematic order 
parameter which is +1 (-1) for a nematic polarized 
along the x (y) direction. Then the symmetry properties 
of the two order parameters show thatO 



sp- 



(26) 



A detailed Landau theory and fluctuation analysis of 
the interplay between these two order parameters is pre- 
sented in Appendix It is our picture, based on the 
physical arguments above, that ^'sp is the primary order 
parameter, and that is tied to it via (p^). 

An important ingredient-iift, the work of Carlson, Or- 
gad, Kivelson, and EmeryllaCj is the crossover from the 
physics of a one-dimensional electron gas moving in the 
vertical direction in Fig |^b at short scales, to a coupled 
two-dimensional system at long scales. This should be 
contrasted from our approach, in which wc do not find 
any quasi-onc-dimensional regime. Although the lattice 
symmetry of our states is similar to those considered 
by these authors, the physics is always intrinsically two- 
dimensional, albeit with a spatial anisotropy. Indeed, the 
arguments for confinement and bond-centered charge or- 
der rely crucially on the two-dimensionality of the sys- 
tem. |— 1 

Zaanen and coUaboratorsO have recently given a bold 
physical picture of the microstructure of stripes and 
their relationship to high temperature superconductiv- 
ity. They take a solitonic, low-dimensional perspective, 
in which special attention is paid to excitations at various 
boundaries and dislocations in the stripe order. In con- 
trast, our perspective here is a higher (two-) dimensional 
point of view, appropriate to systems not too far below 
their upper-critical dimension. So e.g. we view the spin- 
Peierls order parameter "ifsp as a "soft-spin" field with 
large and continuous variations in its local amplitude, 
rather than a field with 4 discrete possible values in dif- 
ferent regions of space separated by identifiable domain 
walls. The advantage of our perspective is that it al- 
lows a complete (in principle) and systematic treatment 
of the coupling of fermionic degrees of freedom to the var- 
ious order gg*ameters in a Landau-theory like scattering 
frameworkOCa. 
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White and ScalapincS have presented results for 
charge order in DMRG studies of doped antiferromag- 
nets. They find that site and bond centered stripes are 
almost degenerate in their energy. They do observe a pe- 
riod 4 bond-centered charge stripe statoD, but have not 
so far seen a statp with the period 2 symmetry of Fig |l|b. 

Castro NetoES has recently presented an interesting 
analysis of a striped superconductor, coupling transverse 
stripe fluctuations between regions of high and low hole 
density. Notably, he also finds pronounced tendency to- 
wards bond-centered charge order with the symmetry of 
Fig^. n 

Lannert, Fisher, and SenthilO have considered Berry 
phase effects in a Zi gauge theory of correlated electron 
systems. Close to half-filling, they find that these Berry 
phases induce bond-centered charge order in both insu- 
lating and superconducting phases, in a manner closely 
analogous to that found in Ref ^ in paramagnetic Mott 
insulators. They also discuss the critical properties of a 
transition between superconducting and insulating states 
at half filling, with both phases containing a background 
charge order with the symmetry of Fig |l|b. They pay 
particular attention to the influence of the gapless nodal 
fermions on the vortices in the superconducting order. 
However, in the presence of the background charge or- 
der at wavevector (tt, 0) , their nodal fermions (which 
are at (±7r/2, ±7r/2)) should immediately acquire a gap, 
and their theory then reduces to a dualized version of 
familiar bosonic theory of the superconductor-insulator 
transition. 
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APPENDIX A: INTERPLAY OF SPIN PEIERLS 
AND NEMATIC ORDERS 

In Section || we discussed the relationship between 
states with spin-Peierls order, with the symmetry of 
Fig |l]b, and the nematic electtepic states discussed by 
Kivelson, Fradkin and EmeryE^I. We noted that the 
spin-Peierls state was characterized by a complex or- 
der parameteiO (which took the values 1, i, —1, 
—i on the four states obtained by rotating Fig |l|b about 
any lattice site) , the nematic was characterized by a real 
Ising order parameter and the symmetry properties 
of these two order parameters implied the relationship 
(p6|) between them. More generally, we can write down 



the following simple effective action for them, keeping all 
low order terms consistent with the underlying symme- 
tries: 



*S'sp-n 



1/T 



dr I d X 
1 



|a.*.p|2 + c2|V,*sp|' + ri|*,p| 



^(9.$„)' + ;^c2(V.<i>„)2 + !J$2 
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(^4 , 

V sp 



sp/ 



W 



(Al) 



where t is imaginary time, and ui > v > 0, U2 > 0, and 
U2{ui — v) > ■Up' required for stability. We expect that the 
appearance of nematic order will enhance the probability 
of spin Peierls order, and so w > 0. We have, for now, 
neglected the couplings of these order to the fcrmionic 
excitations, and will consider the consequences of these 
later in this section. 

First, let us analyze S'sp-n in mean- field theory. The re- 
sults of such an analysis are shown in Fig 14, There are 
three phases: (i) the symmetric phase, where no lattice 
symmetry is broken, (ii) the spin Peierls phase, where 
both the spin Peierls and nematic order parameters are 
non-zero, and (iii), the nematic phase, where only the 
nematic order is non-zero. There can be a second-order 
phase transition between any two of these three phases, 
as indicated in Fig|lJ. The position of these second-order 
boundaries can be determined by a conventional Landau 
theory analysis in powers of the associated order param- 
eter. In this manner, we find that the second-order line 
between the symmetric phase and the spin-Peierls phase 
is at ri = 0, while that between the symmetric phase and 
the nematic phase is at r 2 = 0. Finally, the second-order 
transition from the nematic to the spin-Peierls phase is 
at 



(A2) 



Notice that these three second-order lines all appear to 
meet at the origin r\ = 0, r2 = 0. However, as indi- 
cated in Fig ^4|, this is pre-empted by a line of first order 
transitions close to the origin. The reason for this may 
be seen by the following simple argument. Imagine in- 
tegrating out the $„ fluctuations to derive an effective 
action for the ^'sp: this always induces an effective quar- 
tic term ~ l^spl'* with a coefficient ~ — 5^/|?'2|- So for 
small enough |r2|, the net coefficient of always be- 

comes negative, and this drives the transitions involving 
onset of non-zero (5'sp) first order. A consequence of this 
line of first-order transitions is that there are two tricriti- 
cal points Ti , T2 , and their positions are indicated in the 
caption to Fig |l^. 

Let us now consider the nature of the fluctuations in 
the vicinity of the second-order phase transitions. We 
examine first the transition from the symmetric to the 
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Ti Symmetric 




FIG. 14. Mean-field phase diagram of the action Sap-n 
in ( |A1| ). The symmetric phase has no broken lattice sym- 
metries. The thin (thick) lines represent second (first) or- 
der transitions. The positions of the second order transi- 
tion lines is indicated in (A2) and above. The tricritical 
point Ti is at n — 0, r2 — 4:g^/{ui — v), while the tri- 
critical point T2 is a t the r2 < which is a solution of 
2u2 q'^ + 4qw y/ —r2 /u2+r2(u2(ui —v) —w'^ /2) — 0, and ri given 
by (A2) for this value of r2. At T = 0, the coupling to the 
nodal fermions in a d-wave superconductor drives the traHfii- 
tion from the symmetric to the nematic phase first ordeiEj. 



nematic phase, which has an Ising order parameter. For 
T > 0, these transitions will remain in the universal- 
ity class of the two-dimensional Ising model. However, 
at T = 0, we have to also consider the influence of the 
fermionic excitations. A general analysis of such effects 
has been presented in Ref |7j, and we recall the results 
relevant to our discussion here. The T = ground 
state is either an insulator or a superconductor, and only 
the gapless nodal fermionic excitations can influence the 
nature of the critical properties of a zero temperature 
quantum phase transition. In an anisotropic, d-wave su- 
perconductor, the four nodal points are described by a 
Dirac-like Hamiltonian for 2, four-component, fermionic 
Nambu spinors 5*1, ^'2 (we will follow the notation of 
Ref ^ The nematic order parameter has zero net mo- 
mentum, and so has_a simple linear coupling to these 
fermionic excitationalj: 



(A3) 



where r^'^'^ are Pauli matrices in the Nambu space 
(again, in the notation of Ref [74[). A renormalization 



group analysis of the consequences of (A3) has been car- 
ried out, and it is found that the couplings to the fermions 
drives the T = transition between the symmetric and 
nematic phases flrst order. 

Now let us consider the remaining two second-order 
transitions in Fig for which is the order parame- 
ter. In contrast to ^n, the spin-Peierls order has momen- 
tum (tt, 0) or (0, tt), and so has no simple linear coupling 
to the nodal points, barring the exceptional case in which 
the nodal points are exactly at {tt/2, 7r/2). In general, the 
couplings between and the fermions are of the form 



(^tr^^i + *2r"*2) 



(A4) 



as can be expected from ( pq ) and (A3), and also of the 
form 



(A5) 



Unlike (A3), the couplings in (A4,Af;) are irrelevant, and 
so the fermions do not modify the leading critical proper- 
ties of any transition in which ^^sp is the order parameter. 

For the transition between the nematic and spin- 
Peierls phases, the order parameter is either the real or 
imaginary part of '^sp- the presence of nematic order in 
both phases therefore makes the residual order param- 
eter Ising-like. So this transition is in the universality 
class of three-dimensional (two-dimensional) Ising model 
at r = (T > 0). 

Finally, we consider the fluctuations near the second- 
order transition between the symmetric and the spin- 
Peierls phase. As we have just noted, the fermionic exci- 
tations can be neglected even at T = 0. Further, as ^'sp 
is the primary order parameter, we can integrate out the 
fluctuations at the cost of modifying the couplings in 
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the effective action for the ^fgp alone. The resulting ac- 
tion has the same form as 5, after dropping all the terms 
involving so near the critical point, we need only fol- 
low the flow of the nonlinear couplings ui and v. For the 
T = transition, the standard Wilson-Fisher analysis of 
a phase transition in three spacetime dimension leads to 
the renormalization group flow equations 



-— = ui - 5ui - 9v 
al 

— = w — buiv, 



(A6) 



where we have absorbed a phase-space factor by a rescal- 
ing of the couplings in a standard manner. The only 
stable fixed point of these equations is Ui = 1/5, v = 0. 
At the point v = 0, the effective action for ^fgp has the 
additional symmetry under the global U{1) change in the 
phase of ^'sp. So the T = critical point is in the uni- 
versality class of three dimensional XY model. Related 
considerations can be applied to the T > transition 
between the symmetric and spin-Peierls phases: there, 
the anisotropy associated with v leads to phase tran- 
sitions with continuously varying exponents, associated 
with physics of the Ashkin- Teller modelLa. 

We conclude this section by summarizing the main 
physical implications of Fig |lj. Starting from the sym- 
metric phase, one may break a square lattice symmetry 
either by lowering temperature, or by reducing doping. 
One route to this is a direct second-order transition from 
the symmetric to the spin-Peierls phase; the latter phase 
also has nematic order. The implications of the physical 
arguments we have presented in Section | is that this is 
the preferred route. Alternatively, we can first undergo 
a transition to the nematic phase, and possibly have a 
further symmetry breaking into the spin-Peierls phase. 



APPENDIX B: COMPUTATIONS FOR THE 
TWO-LEG LADDER 

We begin a bond operator theory of doped two-leg lad- 
der by writing down the Hamiltonian for the pure t-J 
model. The effect of Coulomb repulsion will be consid- 
ered in Appendix 

H = Ht+ Hj 

i 



(Bl) 



where i indicates the position of the i-th dimer along the 
ladder direction and («, j) stands for the summation over 
the nearest neighbor, i.e. j = i + 1. 

First, let us take a look at the hopping Hamiltonian 
Ht- Hopping between the sites inside a given dimer is 
given by: 



(B2) 



Hopping between the sites of adjacent dimers is ob- 
tained as follows. 

X \d]hija + -^Sachlj^Sj - -^Sae^^d^ljdtjP^ 



+ ■^d^Sjeachl.^hlj^ + ■^d]s\eabh2ibhija 



1 



£ac£aeO-cbCred4a^jl3lT'2ibhl^d 



(B3) 



where any term containing a sigle t boson operator is 
ignored assuming that there is no magnetic order. Simi- 
larly, 

^ h\ia^2jadid^j + ]^£ab£acs\s.,hiibh\ 



+ -^d^Sjeach\,^h\j^ + -^d]s\eabhiibh2ja 



(B4) 



By combining the above two equations, we get 

AiaClja + 4mC2ja + {h.C.) 
= did]{h\i^hija + /lL^2ja} 

+ ]j^s\^j^abeac{hlibh\j^ + h2ibh\j^) 
+ -^d^js\eab{hlibh2ja + h2ibhija] 

+ ]^£ac£ae(Jcb^^^J[jjp{hlibh\^^ + h2tbhl^^} + {h.C.) (B5) 

Assuming s and d bosons are condensed, i.e. Si = s 
and di = d, one is able to express Ht as follows. 
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V2tds £ab{hli^hl.^ 



+ Ht2h2 (B6) 



where 



+ (h.C.)} (B7) 



Now let us turn our attention to the Hamiltonian Hj. 
Contribution from the coupHng inside the dimer has been 
shown previously. That is, 



3 1 
Sli • 82, = —-S^Si + -ti^tia 



(B8) 



Coupling term between the spins of different dimers is 
given by: 

SliaSlja = SuaSlja + '^(^ab^cd.fAia^iibhljf.hijd, 

+ -jjAiaKb^libSlja + -jh\ja'^2b^^libSua- (B9) 

Remember that Sua = ■^{4^ia + ^la^'^t -ifai3-,t]ptij) and 

S2ia = -5 + tj^Si + icafjjtlf^tij). The cross term 
can be ignored assuming that there is no contribution 
from the term containing a single t operator or t]ytp with 
a ^ (3. In other words, it is assumed that there is no 
magnetic ordering. Therefore, 



= SiictSija + S2iaS2ja + "J^^h* 



(BIO) 



where 



+ h\iJi2ibh\iJi2jd} 



'■2jc"'2jdS (Bll) 

Coupling terms involving S are evaluated below. 

Ija + S2iaS2ja) 



Therefore Hj is summarized as follows. 



+ 2 X]{4a^ja^»*J "I" ^ia^ja'sjs] 

\ ^ t t U" 
^ V ' 
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XJ 

"2 



i?t4 + iJ^4 



(B13) 



where iV is the number of dimers and s boson is con- 
densed so that Si = s. The total Hamiltonian (without 
constraints) is given by the sum of Ht and Hj. 

H = Ht + Hj 

= -tY^{h\i^h2ia + hli^hiia] 

i 

2ja + h\ja^2ia} 

- ^^2tdsY,£ab{h\iAjb + hliAjb 

+h2jbhlia + hijbhzia} 
i 

+ Ht2h2 + H,. + Hh^ (B14) 

As one can see in the following sections, a convenient 
simplification is obtained by expressing the Hamiltonian 
in terms of the bonding and anti-bonding fermionic op- 
erators which are defined as follows: 



v2 



(B15) 



1. Mean-Field Hamiltonian without Quartic Terms 

In this section wc will ignore the terms from j4 , ff ^4 
and Hf2f^2 . Then, the Hamiltonian is given in momentum 
space representation as follows. 
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+ COS fcj, {ti{ky)ti{-ky) + ta{ky)toc{-ky)} 

ky 

+ ^itis"^ - 2(F) cosky - t}hl^{ky)h+a{ky) 
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+ - cos ky + t}hl^{ky)h-aiky) 

ky 

— 2v^t(Jsy~^ COS ky{h\_^{ky)h\_i{—ky) 



-h}_-^{ky)h)_^{-ky) + (/i.e.)} 



(B16) 



As mentioned previously, it is necessary to impose cer- 
tain constraints. In this study, constraints are intro- 
duced in Lagrange multiplier method. In other words, 
"constraint Hamiltonian" [He) is added to the original 
Hamiltonian. And then the unknown parameters are de- 
termined by the saddle point condition of the ground 
state energy. 

i 

+dldi - 1) 

i 

^ -IJ'Y{ti{ky)ta{ky) + + 2X-1) 

ky 

~ {hla{ky)h+a{ky) + hl^{ky)h-a{ky) 

ky 

-|-2(? - 2x) (B17) 

where x is the hole concentration. 

Therefore the total Hamiltonian with constraints is 
given by: 

H= Neo 

+ Y.iA'^y^i{ky)ta{ky) 

ky 

+Bky {tUky)tU-ky) + ta{kyK{-ky))} 
+ Xl{e+(fcl/)/i+a(fc?/)'i+o(fc3/) 

- D{ky){hl^{ky)hl^{-ky) + h+i{-ky)h+-({ky))} 

+ I^{e-(fcy)^-a(^y)/l-a(/Cj,) 

ky 

+ D{ky){hl^{ky)hl^{-ky) + h . ky)h {ky))} (B18) 

where 

3 

eo = --^Js'^ - -(P + 2X-1) 





- £(2(P - 2a;) 


(B19) 


H,y 


/ 1 1 \ T h* 

~ — fJ/ ^ AJ O CUb /iy 


I JDZU ( 


Bky 


= —s cosky 
z 


(B21) 


e+{ky) 


= t{s'^ -2dF) COS ky-t-£^ 


(B22) 


e-{ky) 


= t{s'^ -2dF) cosky + t-^ 


(B23) 


D{ky) 


= 2\/2tds cosky 


(B24) 



The above Hamiltonian can be diagonalized by using 
Bogoliubov transformation. That is, 

7a{ky) = ut{ky)ta{ky) + vt{ky)tl{-ky) (B25) 
P±a{ky) = u±{ky)h±a{ky) + v±{ky)sabhl_^{-ky) (B26) 



where 



Bk 

Ut{ky)Vt{ky) = 

^k„ 



and 



(B27) 
(B28) 
(B29) 
(B30) 

I' +sei©) c*^') 



1 A , e±{ky) 



u±{ky)v±{ky) = T 



Djky) 

Cl±{ky) 



(B33) 

n^{ky) = ^eUky)+D^{ky) (B34) 

In terms of Bogoliubov variables, Hamiltonian is writ- 
ten as follows. 

H=Neo + J2'^kyli{ky)7a{ky) + ^Yl ('^''v ~ "^^y) 

ky ky 

+ Y^+iky)PUky)Mky) 

ky 

+ Y [^+{ky) - ^Jel{ky)+D^{ky)^ 

ky 

+ Y^-{ky)Pla{ky)P-a{ky) 

ky 

+ Y\^e.{ky)-^el{ky) + D^ky)'^ (B35) 
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a. Ground State Energy and Saddle-Point Equations 
Now the ground state energy per particle is given by: 



<^gr 



N 



Jo 

/"^ rik 

Jo 



I 3 dk. 



= e, 



2 ./o TT 



- r^n^iky)^ r^n.iky) (B36) 

Jo ^ Jo 



where 



Co = —Js'^ - fJ-is"^ -d^ + 2x-^)-^{2 + 2cP - 2x). 

(B37) 

Note that constants obtained from integrating Ai-y and 
e±{ky) are ignored. 

The term contributed from the triplet bosons (ta) can 
be intergrated explicitly by using the elliptic function E. 



_ r '^^n+{ky)- r ^n_{ky) (B38) 

Jo Jo 



where 



4 J 



(B39) 



Parameters in the Hamiltonian are determined by the 
saddle point condition in the mean-field approximation. 
In other words, 



dCgf dCgr d^gj- d^gj- 



(B40) 



ds dd dS, 

Explicit saddle-point equations are written as follows. 



1 __j,f I 2. 



vT+~n VV 1 + ?7 



t/J ['^ dky j st+{ky) + ^/2dD{ky) 

COS ntqi \ 



-se-{ky) + V2dD{ky) \ ^ 

n.{ky) { 



(B41) 



9(/i/J) 2 



(B42) 



dd^ J Jt 



t/J r dky j 2de+{ky) - V2sD{ky) 



^2de.iky)-V2,Diky)\^ (B43) 

U-{ky) I 



dfcj, / c+{ky) e-{ky) 



TT if^+(A;j;) ^-{ky) 



(B44) 



6. Limiting Case with d = Q 

The above saddle-point equations are highly non-linear 
as a function of parameters s, ^i, d, and ^. Since analytic 
solution is not accessible in this case, it is natural to make 
use of a numerical method. In general, solving a set of 
non-linear equations numerically is quite sensitive to the 
initial guesses of parameters. Therefore, it is very desir- 
able to have some limiting situation where an explicit, 
analytic solution is available. Such a situation is realized 
when we set d = 0, i.e. turning off the pairing terms. 
Then, the ground state energy is given as follows. 



dk 

/ ^{e+{ky)-n+{ky)} (B45) 
Jo 



+ 



where the contribution from the anti-bonding fermion 

is ignored because, in the case of t/J > 1, the 
anti-bonding fermion will have a reasonably larger en- 
ergy than the bonding one, separated by roughly 2t. And 
therefore its contribution to the ground state energy is 
very small. 

The fermion contribution can be simplified further. 



rdk^ 
J 



^L^[ky)-^+{ky)] 

2 r 

= - / b+{kv)\d=o 
2 

= — I dky {ts^ COS ky — t — ^) 

Jkp 



-- {{t + Oi'^ - kF)+ts'^smkF} 
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— —ts^ {ttx cos (ttx) — sin (ttx)} (B46) 

TT 

where the chemical potential ^ and the Fermi wavevector 
kp are determined by: 

^ h\^{ky)h+a{ky) = X ^ ^(tt - fcf) = 2x (B47) 

ky 

e+ikp) =0^^^ts'^ COS kp-t. (B48) 

Now the saddle-point equations are reduced to the fol- 
lowing equations. 



1 



3 u 

4 + 7 

=K 



3A 



2^ 
1 + 77 



VTT^ V V 1 + ?? 
— (ttx cos (ttx) — sin (ttx) ) = 

TT 

2 27r I ^ ' 



(B49) 



2?7 



1 



=K 



2r/ 



1 + ?/, 
(B50) 



Note that the other saddle-point equations are explicitly 
solved. In other words, 



d = 

^ = —ts^ cos [ttx) — t 



(B51) 
(B52) 



C onve nient simp lification is realized from the fact that 
Eq.(B49) and (|B5C| ) can be combined to reduce to a equa- 
tion containing only one parameter 77. That is, 



1 (ttx cos (tto;) — sin (ttx)) 

2A TT J 



^-2.-^^K 
2 TT VTT^ 



2ri 
I + 11 



(B53) 



where, once again, 77 — y/T^^JT/I' 

2. Mean-Field Hamiltonian with Quartic Terms 

In this section, effects of quartic terms such as 77^4, 
Hj^i, and -^42 ^2 are included approximately by using 
quadratic decoupling. Let us start from H^i . 

XJ 



XJ 



Y.^iLU(4o.hP - tlUpthtjc} (B54) 



Quadratic decoupling is carried out as follows. 

- {t\jio.){t]pUp) (B56) 
Therefore, Hti is given by: 



— — cos ky {2Pyt\{ky)ta[ky) 

ky 

- Qy {ti{ky)tl{-ky) + U-ky)t^{ky)) } 



(B57) 



where 



Py = (4tj„) = ^ COS ky{ti{ky)t^iky)) (B58) 

ky 

Qy ^ (44) ^T.'^''^''y(t^Mtli-ky)). (B59) 



Note that the translational symmetry is assumed. Now 
let us turn to the Hamiltonian Hf2f^2. 



+ [h.c.)} 



(B60) 



Remembering that, in the present approximation, non- 
zero quadratic contraction is obtained only when a = (3 
for {tl^tp), and a = b for {h^^ht), the Hamiltonian iJt2^2 
is quadratically decoupled as follows. 

+ EacSaeCr^b^Hd 4jjf3{hlibh\j^ + {h.C.)} 

~ eac£aeCr"b^eb Ajj-^ihubhl^f, + {h.C.)} 

^{tU3o){h\^ahua) + {h.C.) 
= -Py (hlja^T-lia + /iLa^lja) ^ (tl^tja + ^L^ja) 

+2PyIly (B61) 
where 

Uy = {h\^,hua) =Y''°^hKa{ky)hia{ky)). (B62) 
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Once again it is convenient to use the bonding and anti- 
bonding fermion representation. By using the reflection 
symmetry, 11,^ is given by: 

fey 

= Y, COS ky{hl^{ky)h2a{ky)) 

ky 

= COS ky{{hl„{ky)h+a{ky)) 

ky 

+{hlaiky)h-a{ky))} (B63) 
Finally, the Hamiltonian H^2f^2 is given by: 

Ht2h2 = tPyY^cosky (hl^{k)h+a{k) + hl^{k)h-a{k)) 

+ 2ttly^^ cos kyt^„{ky)ta{ky) 



- AtNPylly 



(B64) 



Now, quadratic decoupling of If/j4 is performed. First, 
Hh4 is given by: 



ijd 



(B65) 



Suppressing the site index within a dimer, i.e. 1 or 2, let 
us first consider the following term inside summation. 

3 3 
= -^^y{h]ahia + hlhja) + -^l 

+6A2 (B66) 
where, once again, IIj, is given by: 

= \Y^cosky\^h\SK)h+a{ky)) 

ky 

+{hla{K)h-a{ky))y (B67) 
and also Ay is defined as follows: 



(B68) 



+{hU{ky)hl^{-ky)) 



Therefore H^i can be written as follows. 

Hh^ = -^AJ^cosfc3^|nj^/i^,_„(fcj^)/i+a(fcy) 

k ^ 

rvy 

+ 2Ay (hl^{ky)hl^{-ky) + h + ^{-ky)h+^{ky] 

- \>^J ^cosky[llyh}_^{ky)h-a{ky) 



+ 2Ay (hl^{ky)hl^{-ky) + h - ^{-ky)h {ky) 



+ -7VAJ(n^+4A^) 



(B69) 



Including all the quartic terms, the full, mean-field 
Hamiltonian is finally given by: 

H=Nea + Y.{'^k/Mta{ky) 

k ^ 

+ B'ky {tUky)ti{-ky)+t4ky)t4-ky)) ^ 

+ Yl \^^'+iky)hlaiky)h+a{ky) 

ky 

-D'+{ky) (^hl^{ky)hl^{-ky) + h+i{-ky)h+i{ky) 
+ J2\^'-^''y^f'-aiky)h-a{ky) 

k ^ 

fvy 

+ D'_{ky) (^hl^{ky)hl^{-ky)+h.l{-ky)h.^{ky)) ^ 

(B70) 



where 



3 



eo = --^-J^^ - - cP + 2a; - 1) 
-^(2fP - 2x) 

Af.^ = Ak^ + XJPy cos ky + 2tl{y COS ky 

= ^ — IJ, + XJ{s^ + Py) COS ky + 2tlly COS k 

B'ky =Bky - ^Qy COS ky 

= -^i^^ - Qy) COS ky 

I 3 

^+{ky) = ^+{ky) — -AJIlj^ cos ky + tPy COS ky 

3 

= t{s'^ — 2(? + Py) COS ky — -A Jlly COS ky —t— 

I 3 

£-{ky) = e-{ky) — -AJIIy COS ky + tPy cos ky 

3 

= t(S^ — 2(? + Py) COS ky — ^A Jlly COS ky + t — ^ 

' 3 

D^{ky) = D{ky) + -XJAyCOSky 
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= 2\f2tds cos ky + - A JAy cos ky 

f 3 

D_{ky) — D{ky) — -XjAy COS ky 

3 

= 2\/2tds COS ky — -XJAy cos ky 



(B71) 



Note that all the constants which are not explicitly de- 
pendent on s, fi, d, and ^ arc ignored. As previously, the 
above Hamiltonian can be diagonalized by using Bogoli- 
ubov transformation. That is, 

H = Neo + Y.^k,'yi{kyha{ky) + ^ Y.i'^'ky - ^'ky) 

+ Y.^'+{ky)Pia{ky)P+a{ky) 

ky 

ky 

+ Y,^'_{ky)PUky)P-a{ky) 

ky 

+ Y,(eUky)-^lUky)) 

ky 

(B72) 

Saddle-point equations of the full Hamiltonian with 
quartic terms are similar to those without quartic terms 
except the fact that even the t boson part is not expressed 
explicitly. 



4 J 2 Jo TT 

■ COS ky 



t/J dky ( se'^iky) + V^dD'^{ky) 







^+{ky) 



^ se'_{ky) + V2dD'_{ky) , _Q 



n'_iky) 

d{egr/J) M 



3 r dky^ 
2 io 7^ '^fc. 



(B73) 



(B74) 



^ -2-- 



J "Jt 

+ -U^ / COS k 



t/J r dky ( 2de'+{ky) - V2sD'^{ky) 

~ / COS rC-ji \ ~~, — ' ^ 

io n ^+{ky) 



2de_{ky)-V2sD'_{ky) \ _ 
Q'_{ky) ^ 



(B75) 



dky J e_|_(fcy) i_{ky) 



n dky I e+(ky, 

Jo n \ 0.'_^{ky 



+ I —< ^t;.°; + ^7;/; > = o (b76) 

However, the above saddle-point equations are not 



j) ^-{ky) 
Ut 

meaningful without the knowledge of Py, Qy, Hy, and 



Aj, used in A^.^, Bj,^, e^{ky), and D^{ky). These are the 
formulas for them. 



3 r dky 



COS ky —r" 







Qy = -'i I ^COSky-p- 



f dky 
1 r dk„ 



(B77) 
(B78) 



e+(fcj;) i-{ky) 



Uy = -ll tcos..|^.gg^^ (B79) 



A -i r^^cosfc 

2I0 TT ''°'''M 



D'+{ky) D'_{ky) 
{ky) Vl_{ky) 



(B80) 



Upon solving the saddle-point equations, it is conve- 
nient to view Py, Qy, Hy, and A^ as unknown parame- 
ters similar to s, d, /z, and ^. Then, the original saddle- 
point conditions are equivalent to solving four saddle- 
point equations and four equations for Py, Qy, Hy, and 
Ay simultaneously. 



APPENDIX C: COMPUTATIONS FOR THE 
SQUARE LATTICE 

This section is devoted to the two-dimensional array 
of coupled t- J ladders. Hamiltonian for the coupling be- 
tween two-leg t-J ladders can be written as follows. 

i {n,ra) 



i (n,m) 



where the index i indicates the i-th dimer within a given 
ladder and the indices n {m = n + 1) indicate the n-th 
(m-th) ladder. Derivation of the mean-field Hamiltonian 
for the dynamics perpenticular to the ladder direction, 
is rather similar to that of the parallel direction. But 
caution should be used in keeping track of signs of various 
terms in Hamiltonian. Let us start from the hopping 
term between ladders. 

ft If 1 a t \ 



1 



+ -^dsSachlinahlimc + -^dsE abhunbhu 
7^Sac^aeO'f^0^jt\^^timphiinbh\^^^ 



(C2) 
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Also, assuming that there is no magnetic ordering, the 
product of spin operators is given by: 

= SiimaS2ina + '^^ab^cd^lima^Umbh^inc^^ind 

=^ (tlma^ina + t\na^ima + A.m<xA,noi timatina) 

+ ^'^ab'^cd^lima^limbhli,ich2ind (C3) 

Performing the quadratic decoupling and ignoring ir- 
relevent constants, the Hamiltonian H± is writtten as 
follows. (Note that is the momentum associated with 
the perpendicular direction to ladders, while ky is with 
the parallel direction. Also, in laboratory unit, momen- 
tum in the parallel direction is given by py = ky/a, and 
that in the perpendicular direction is given by p^; = kx/2a 
where a is the lattice spacing and —wKkx, ky < n.) 

H± = HJ_f^2 + Hj_t2fi2 + Hj_t2 + Hj_t4 + HJ_f^4 (C4) 



H^,h- = -{s^ - 2(P)^cosA;4/.^„(k)/i+<,(k) 
k 

-/iL„(k)/i_„(k)} 

- V2t ds^coskx{h\_^{\s.)h\_^{-k) 

k 

- V2t ds^^cosk^lhl^Ck^hl^i-k.) 

k 

+/._i(-k)/i_T(k)} (C5) 
H±,t^h^ = --P^^coskx{h\_^{k)h+a{^) 

k 

-/iL,(k)/^_„(k)} 

- i ^ coskxtl{k)ta{k) 



AJ. 



H±,t- = -^s2^cos/c424(k)f„(k) 



+4(k)tt(-k)+f„(k)t„(-k)} 



(C6) 



(C7) 



A J 



if±,t4 = — ^cosfc,{2P,4(k)t„(k) 

k 

-Ox {ti{k)tl{-\c) + t„(k)t„(-k)) } (C8) 



H±^h* = --A Jn^ Vlcosfca;{/i^(k)/i+a(k) 
^ k 

-/iL,(k)/i_,(k)} 

- -A'jA^^cosfc^{/i^^(k)/i^^(-k) 

k 

+/l+i(-k)/i+T(k)} 

+ -A JA^^cosfc^{/il^(k)/iL|(-k) 

k 

+/i_i(-k)/i_T(k)} (C9) 

where 

f'x = ^COSfc^(tt(k)i„(k)) 
k 

Q, = 5]cosfc,(tt(k)4(-k)). 

k 

= i^cosfc, {(/it„(k)/i+„(k)) - (/il„(k)/i_„(k))} 

A, = i^cosfc,|(/i;^(k)/i^^(-k)) 
-(/iLT(k)/il^(-k))| 

(CIO) 

where the bonding and anti-bonding fermions are as- 
sumed to be well separated in energy so that the mixing 
terms between them are ignored. 

Combining H± with the Hamiltonian describing the 
dynamics inside ladder, the mean- field Hamiltonian for 
the two-dimensional array of coupled t-J ladders is ob- 
tained. 

H = Nho/2 + ^ {ik4(k)ta(k) 

k 

+Bk {tUk)ti{-k) + t„(k)t„(-k)) } 

+ 5^{e+(k)/i^(k)/i+„(k) 

k 

-D+{k){hl^{k)hl^{^k) + h+i{-k)h+^{k))} 
+ ^{e_(k)/iL„(k)ft_„(k) 

k 

+D_{k){hl^{k)hl^{-k) + /i_i(-k)/i_|(k))} (Cll) 
where 

3 



eo = -^Js^ - m(s^ -cP + 2x-1) 
-C(2(P - 2x) 

Ak = -J — /X + A Js^ ( cos ky — cos kx 
4 \ 2A 



+ XJ\Py COS ky -j- COS kx 



(C12) 
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+ 2t I Ily COS ky — -2^^x COS k. 



By. = S COS ky COS kx 

2 \ ^ 2A 



(C13) 



^ [QyCosky + ^QxCoskx\ (C14) 



A 



e+(k) = t(s^ — 2(P) I COS ^v'^ ^x] —t — ^ 



5a. n.c„,*.4n,cos.; 



t I Py COS ky — —Px COS fc 



2t 



(C15) 



e_(k) = t{s^ — 2(P) [ cosky — — cos fc^; j + t — 



3 / 

-A J I n„ cos fc„ — —Ha; COS 



.y uuo fby ^^±±x cob ft-x 



+ i I Pj/ cos fcj; + —Pa; COS kx 



(C16) 



£>+(k) = 2\/2tds l^cosfcj, + — cos/c^J 

+ ^ A J I Ay cos fcy + ^ Ax cos A;^ j (C17) 



^_(k) = 2\/2tds l^cosfcy — ^ cosfcxj 

3 / a' \ 

- -A J I Ay cos ky - — Ax cos kx \ (C18) 

Diagonalization via Bogoliubov transformation and 
the correspondent saddle-point equations are similar to 
those of single ladder case. That is to say, 



-^w(k) 



-/ /5^Mk)-/ /5^Mk) (C19) 



where 



^(2 + 2(P - 2x) 



(C20) 
(C21) 



fi±(k) = ^el(k) + £>|(k). (C22) 
where the upper and lower limit of integral, tt and 0, are 



not explicitly written. Saddle-point equations are given 
as follows. 

d{egr/J) _ _3 _ M 
~ 4 J 



'2 J J -^^ 

ik - 2Bk 



X 



t/J f f <P\^ ( , t , ^ 
-—] J -^\oosky + -coskx 



X 



se+(k) +72rf£'+(k) 
^//5fcos.y-^cos.. 



_ 3 /■ /rf^A ^ 
2 y 7 TT^ (Dk 



9(^2 J J i 



-T cosky + -coskx 



d J J TT^ « ' 2t 

2de+{\s.) - y2sL>+(k) 
^^ 2de-(k)-V2g-i^_(k) _^ 



(C23) 



se_(k) + y2dD_(k) „ 
X = — = 

fi_(k) 

■^'"^^•^ " (C24) 



(C25) 



fi_(k) 

9('9r/J) _ 2 2^-2 2^ 

^ ^ — I ^ (C26) 



lf]+(k) 17- (k). 

Also, the parameters due to quartic terms are defined as 
follows. 



'~2j J TT^ '^°'''^a;k 



' d^k , Bk 

1 r fd^k 



Qy ^ I I cosky — 



(C27) 

(C28) 



n„ = -- / / ^ cos/c^ { ' + , ^ ' \ (C29) 



1 f f d^k 



' j+(k) ^_(k) 

n+{k) fi_(k) 



(C30) 
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cos fea; — 



TT' 

7r2 



COS Kx 

Wk 



7r2 ll7+(k) f]_(k) 



COS kx 



D+{k) ^ ^_(k) 
f2+(k) f2„(k) 



(C31) 
(C32) 

(C33) 

(C34) 



As previously, it is convenient to have some limiting 
case when the explicit solution is available. We will take 
the case of a; = 0, i.e. no-doping case. 



3 M 3 
"4~ 7+2^ 



d{egr/J) 



djegr/J) _ 5 _ _2 _ 3 
d{n/J) ~ 2 ^ 2 



rf'k /k 



VI + Vh 

d^k 1 + f/k _ 







(C35) 
(C36) 



where 7] = i/l^fj^/ 7 and /k = cos fcy — ^ cos fcy . As 
before, it is possible to combine the above two equations 
to get a single equation only for r/. That is to say. 



ry 5 3 
2A " 2 ~ 2 



//^ 



1 



(C37) 



1. Superconducting Order Parameters 

This subsection relates the parameters in the above 
computations to the physical order parameters describing 
the nature of the superconducting state. In the case of 
two-dimensional array of coupled ladders, there are three 
order parameters to be defined: 



XO = ^ab'\Ciina'^2inb) 
X\\ = ^ ab{(^i^iJ^i^naA^inb) 
X-L = ^ab{c[i^n+l,a^2inb> 



(C38) 

(C39) 
(C40) 



+ -j^eac£be(Tcd'^e}{^l+l,naAn0^'^,i+^^ndh2inf)} 

= {(F- + Q)r - V2dsii 



(C42) 



where Q and F are defined previously, and 11 is given by: 

e-{k,q) 



n = -l / /^cosJi±(^- 

2j J 7r2 \n+{k,q) 



n_{k,q) 



(C43) 



Similarly, pairing order parameter in the perpendicular 
direction, x±, is computed to give the following formula. 



X± = £ab{(P {hl^naf^\i^n+l,b) + l^^ ac£bdS^ {hunch2i,n+l,d) 
+ -^£bcds(4i„a^2i,n+l,c) + ■^eacds(/Hi„c/lH,„+i,b) 

- ~eace6eO-cdCrf^(4na4,n+l,/3^1»nd/i2i,n+l,/)} 



(C44) 



where Q±_ and r_L are defined previously, and Hj, is given 
by: 



- If [ dkdq ( e+{k,q) e-{k,q) 

n_L = - - / / — TT- cos q i 1- 

2j J 7r2 \n+{k,q) n_{k,q) 



(C45) 



APPENDIX D: EFFECT OF COULOMB 
INTERACTION 

This section is devoted to the effect of Coulomb re- 
pulsion. Hamiltonian due to the Coulomb interaction 
between the nearest neighbors is written in terms of the 
bond operators as follows. 



The first one, xoi is the pairing order parameter within 
dimer. Using constraints, it is written as follows. 

1 1 

XO = £ab{{h\ina(^in + -^acs\jX2inc - -^£ac(^cd4naf^'^ind} 



Hcoulomb = V ^ ^ nii 



<,n2v. 



{^linb^i"^ + ^^bcslnhiinc + ^^bcf^cd^lnffhlind}) 

(C41) where 



+ ^ ^{«lm"ljn + n2inn2jn} 

(.id) " 

+ ^ ^ n2innum 



(Dl) 



= £ab ■ -^(£acsLrfm/l2mc4m6) = "^^^ 



The pairing order parameter in the parallel direction, 
X\\, is evaluated as follows. 

X|| = Sab{(P {hl^i+i^nal^linb) + ^^acEbds"^ {h2,i+l,nah2ind) 



i {n,m) 



= 1 — h2ina^2ina — din'^in 



(D2) 



and 
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in — '^2ina 



= - h\^^^hiina - d\„din. (D3) 

The first term is evaluated as follows. 

i n 
i n 

i n 
i 



I n 
t2 



-{l-2x + d^) 



(D4) 



where N"^ /2 is the total number of dimers. Note that the 
constraint conditions are used. 

The second term is due to the Coulomb interaction 
between the dimers along the ladder direction, which is 
computed by: 

{id) " 
(ij) " 

= N^{l-(P){l-2x + (P) 



Similarly, the third term is evaluated as follows. 



(D5) 



(D6) 



Therefore, 



iV2 



^E E ^\ina^linah\i^^h2irab 



i (n,m) 

— y(4 - 3(P)(1 - 2a; + <?) + ;,4 



(D7) 



The full treatment of HqouI, h'^ is beyond the scope of 
our present paper. However, we will follow the philoso- 
phy of BCS theory in that electrons interacting with the 
Coulomb repulsion under the neutralizing background 
are described in terms of the Landau-Fermi liquid. As- 
suming that the hopping parameter, <, and the spin- 
coupling parameter, J, are already renoramilzed quati- 
ties due to the Coulomb repulsion, one has only to con- 
sider the Hartree contribution of Hcoui /i* i which is given 
by: 



(D8) 



which is obtained by using the constraint conditions. 
Also, note that the mixing between the bonding and anti- 
bonding fermions is ignored. 

Therefore, in our approximation, the Hamiltonian (per 
site) due to the Coulomb repulsion is given as follows: 



He 



oulomb 



V{A - 3(F){1 - 2x + (P) + Hl 



= V{A^8x + 3x'^) + V(F (D9) 

Since Eq.(p9|) depends only on cP, the modification of 
saddle-point equations is straightforward. That is. 



dd^ 



= 



d{egr/J) 



V = Q 



(DIO) 



All the other saddle-point equations are the same as be- 
fore. 

Now it is reasonable to assume the limit of large 
Coulomb repulsion which not only captures the essential 
physics that we are interested in, but also will simplify 
explicit computations. So we will take the limit V — > oo; 
consequently c? = 0. While this limit is not physical for 
large doping, our primary interest lies in the case of small 
doping. 
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